PRECIPITACIÓN EN PUEBLA

Precipitación v4:
https://www.ncei.noaa.gov/products/land-based-station/global-historical-climatology-network-monthly https://www.ncei.noaa.gov/data/ghcnm/v4/precipitation/ https://www.ncei.noaa.gov/pub/data/ghcn/v4/products/StationPlots/MX/

Nombres:
https://www.ncei.noaa.gov/pub/data/ghcn/v4/
https://www.ncei.noaa.gov/data/ghcnm/v4/precipitation/doc/

Ubicación de la estación: 19.0000, -98.1833

CONFIGURACIÓN DE LA NOTETBOOK¶

In [27]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
In [28]:
sns.set(style="whitegrid")
In [29]:
# The color palette is made up of the 20 colors. Hex color codes:  #395e77, #413344,  #614c65,  #806485,  #936397,  #a662a8,  #664972,  #463c57,  #6e8da9,  #91bcdd,  #567d99,  #305662,  #264d4d,  #315c45,  #8a9a65,  #b6b975,  #b65d54,  #b60033,  #98062d and  #800022.
# https://colorkit.co/palette/413344-614c65-806485-936397-a662a8-664972-463c57-6e8da9-91bcdd-567d99-395e77-305662-264d4d-315c45-8a9a65-b6b975-b65d54-b60033-98062d-800022/
sns.set_palette(["#194a7a", "#b60033", "#315c45", "#b65d54", "#B9445F", "#567d99", "#395e77", "#413344", "#614c65", "#806485", "#936397", "#a662a8", "#664972", "#463c57", "#6e8da9", "#91bcdd", "#305662", "#264d4d",  "#8a9a65", "#b6b975", "#98062d", "#800022"])
sns.color_palette()
Out[29]:
In [30]:
# Personalización global con matplotlib
plt.rcParams.update({
    'axes.titlesize': 14,        # Tamaño del título
    'axes.titleweight': 'bold',  # Negrita en el título
    'xtick.labelsize': 8,        # Tamaño de los xticks
    'ytick.labelsize': 8,         # Tamaño de los yticks
    'grid.color': 'gray',         # Color de las líneas del grid
    'grid.linestyle': '--',       # Estilo de línea (puede ser '-', '--', '-.', ':')
    'grid.linewidth': 0.5,        # Grosor del grid
    'axes.grid': True,            # Asegura que el grid esté activado
    'axes.grid.axis': 'both',     # Aplica el grid a ambos ejes
    'lines.linewidth': 1.2,       # Grosor de las líneas
    'figure.figsize': (12, 6),   # Tamaño de la figura
})

CARGAR DATOS¶

In [31]:
data=pd.read_csv('MXN00021035.csv')
In [32]:
pre=data.iloc[:,6]  # Precipitacion, es la columna 5
date=data.iloc[:,5] # Date, es la columna 6
date = date.astype(str).str.replace(r'(\d{4})(\d{2})', r'\1/\2', regex=True)    # La fecha está como 195210 y la pasamos a 1952/10 
date = pd.to_datetime(date, format='%Y/%m')                                     # Lo convertimos en fecha
pre = pd.Series(pre.values, index=date)                                         # Creamos una Serie
pre
Out[32]:
195209
1952-10-01     174
1952-11-01     592
1952-12-01       0
1953-04-01      47
1953-05-01     137
              ... 
2009-08-01    1474
2009-09-01    3161
2009-10-01    1190
2009-11-01     155
2009-12-01      51
Length: 672, dtype: int64
In [33]:
import matplotlib.dates as mdates
In [34]:
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3))  # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha
plt.plot(pre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Precipitación Mensual en Puebla")
plt.savefig('imagenes/01-01-precipitacion.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [35]:
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12))  # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha
plt.plot(pre)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Precipitación Mensual en Puebla (2000-2010)")
plt.xlim(pd.Timestamp('1990-01-01'), pd.Timestamp('2010-01-01')) 
plt.savefig('imagenes/01-02-precipitacion-zoom.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

TRAIN Y TEST¶

In [36]:
# Partir la serie para train y test
pre_total = pre.copy()          # Copia de la serie original

# Todas hasta los ultimos 12 meses
pre = pre_total[:-12]           # Entrenamiento: todos menos los últimos 12 meses
pre_test = pre_total[-12:]      # Test: últimos 12 meses

SERIE SIN TRANSFORMACIONES, SOLO ESTANDIARIZADA¶

Como la vamos a comparar con la transformacion de Yeo-Johnson, estandarizamos para comparar los modelos.

In [37]:
mu = pre.mean()                # Media
sigma = pre.std()              # Desviación estándar

# Normalizar la serie
zpre = (pre - mu) / sigma    # Estandarizar la serie

ESTACIONARIEDAD¶

In [38]:
from statsmodels.tsa.stattools import adfuller  
In [39]:
adfuller(zpre)
Out[39]:
(-6.150368236022414,
 7.588129761677003e-08,
 11,
 648,
 {'1%': -3.4404817800778034,
  '5%': -2.866010569916275,
  '10%': -2.569150763698369},
 1333.5800521506967)
Estadístico ADF         = -6.150  
Valor-p                 = 7.59e-08  
Número de rezagos       = 11  
Número de observaciones = 648  
Valores críticos:
    1%  -> -3.4405
    5%  -> -2.8660
    10% -> -2.5692
Log-likelihood          = 9999.935
  • El estadístico ADF es $-6.150$, menor que los valores críticos a los niveles del 1%.
  • El valor-p es $7.59 × 10^{-8}$ mucho menor al nivel de significancia $\alpha$ = 0.05.

Por lo tanto, se rechaza la hipótesis nula, es decir, se concluye que la serie es estacionaria.

AUTOCORRELACIONES¶

In [40]:
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf, acf, pacf
import matplotlib.pyplot as plt
from fac_y_facps_significativas import *    # Este es el archivo que se usa para el curso
In [41]:
fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(zpre, lags=48, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/02-01-fac-facp-serie.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [42]:
fac = FAC(len(zpre), acf(zpre, nlags=48)[1:] )
Valores de autocorrelacion significativos:
r1: 0.5114347926277153
r2: 0.23788150833700863
r4: -0.31136216826642676
r5: -0.5289967604156934
r6: -0.5765988530134252
r7: -0.4540298620448335
r8: -0.24420551667402526
r10: 0.2960004737880804
r11: 0.5240658344251691
r12: 0.6173499572874946
r13: 0.424401244621075
r14: 0.18265000690415573
r16: -0.31903920298321264
r17: -0.514286815393446
r18: -0.527829591696689
r19: -0.4187324328138189
r22: 0.30067304144252044
r23: 0.5028304240974979
r24: 0.5626631891391944
r25: 0.3994834359095428
r28: -0.34516091598549786
r29: -0.48547226384609826
r30: -0.5033566062877032
r31: -0.38761684940293784
r34: 0.3085645527226526
r35: 0.4965812730150526
r36: 0.553098035132298
r37: 0.36241966295184613
r40: -0.3454697704347722
r41: -0.4674618422218585
r42: -0.46761339396960605
r43: -0.34243418620086286
r46: 0.3414843279972229
r47: 0.48860016886221763
r48: 0.5000880408130988
In [43]:
facp = FACP(len(zpre), pacf(zpre, nlags=48)[1:])
Valores de autocorrelacion parcial significativos:
rho 1: 0.5122108697030228
rho 3: -0.1371024981440388
rho 4: -0.35568794886746186
rho 5: -0.3390668723583037
rho 6: -0.26990992305506756
rho 7: -0.13651941979413404
rho 8: -0.082068919602615
rho 11: 0.19446853660047775
rho 12: 0.25021094781029624
rho 19: -0.0924235210146644
rho 24: 0.12444505947195694
rho 32: -0.08477112361318034
rho 36: 0.11301403546356209

Aunque la serie es estacionaria, dado que las autocorrelaciones simples decrecen muy muy lento, hacemos una primera diferencia estacional.

DIFERENCIA ESTACIONAL¶

La diferencia estacional de orden $D$, con una periodicidad estacional $s$, se define de forma recursiva aplicando el operador $\nabla_s$ múltiples veces:

$$ \nabla_s^D X_t = (1 - B^s)^D X_t $$

donde:

  • $D$ es el orden de la diferencia estacional,
  • $B$ es el operador de rezago tal que $B^k X_t = y_{X-k}$,
  • $\nabla_s^D$ representa aplicar $D$ veces la diferencia estacional.

Para esta serie, el orden de la diferencia estacional es $D = 1$ y la periodicidad estacional es $s = 12$:

$$ \nabla_{12}^1 \text{T} (X_t) = \text{T} (X_t) - \text{T} (X_{t-12}) $$

Sea entonces $Wt = \nabla_{12}^1 \text{T} (X_t) = \text{T} (X_t) - \text{T} (X_{t-12})$

PRIMERA DIFERENCIA ESTACIONAL¶

In [44]:
# Primera diferencia estacional con periodicidad 12
dzpre = zpre.diff(12).dropna() 
In [45]:
# Visualizamos la serie
dzpre
Out[45]:
195209
1954-02-01   -0.044262
1954-03-01   -0.621941
1954-04-01    0.357503
1954-05-01    1.805672
1954-06-01    3.170992
                ...   
2008-08-01    0.762672
2008-09-01    0.396090
2008-10-01    0.048802
2008-11-01    0.039723
2008-12-01    0.000000
Length: 648, dtype: float64
In [46]:
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3))  # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha
plt.plot(dzpre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Primera Diferencia Estacional, Estandarizada, Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-01-primera-diferencia-estacional-precipitacion.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [47]:
adfuller(dzpre)
Out[47]:
(-10.871408688779647,
 1.3655020494523785e-19,
 11,
 636,
 {'1%': -3.4406737255613256,
  '5%': -2.866095119842903,
  '10%': -2.5691958123689727},
 1448.5800337609196)
Estadístico ADF         = -10.871  
Valor-p                 = 1.37e-19  
Número de rezagos       = 11  
Número de observaciones = 636  
Valores críticos:
    1%  -> -3.4407
    5%  -> -2.8661
    10% -> -2.5692
Log-likelihood          = 9952.186
  • El estadístico ADF es $-10.871$, menor que los valores críticos a los niveles del 1%.
  • El valor-p es $1.37 × 10^{-19}$ mucho menor al nivel de significancia $\alpha$ = 0.05.

Por lo tanto, se rechaza la hipótesis nula, es decir, se concluye que la serie es estacionaria.

In [48]:
fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(dzpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(dzpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=1, Estandarizada de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-02-fac-facp-serie-D1.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

Se sigue notando la estacionalidad de periodicidad 12.

In [49]:
fac = FAC(len(dzpre), acf(dzpre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos:
r1: 0.0849542779953831
r5: -0.09455404757011154
r12: -0.414183329502961
r69: 0.13811523521698316
r73: -0.1125610007638793
r82: -0.09936326169869979
In [50]:
facp = FACP(len(dzpre), pacf(dzpre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos:
rho 1: 0.08508558290727707
rho 5: -0.08389897753926695
rho 12: -0.435832204242704
rho 17: -0.12023508334631498
rho 24: -0.30097443805877616
rho 29: -0.10198920198664953
rho 36: -0.15950168604080855
rho 48: -0.14036369491687178
rho 57: -0.0916871776051421
rho 60: -0.12396066058923692
rho 61: 0.13800141862943816
rho 69: 0.10269095820413293
rho 72: -0.09462958643919762
rho 82: -0.08002594037040592
rho 85: -0.0819832551248085
rho 93: 0.0922037887338383
rho 109: 0.09582379018763905
rho 113: -0.09105599732819779
rho 119: 0.0917177196667104
rho 120: -0.1355872480803115

SEGUNDA DIFERENCIA ESTACIONAL¶

In [51]:
# 2 diferencia estacional con periodicidad 12
d2zpre = dzpre.diff(12).dropna() 
In [52]:
fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(d2zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(d2zpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')


plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=2 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-03-fac-facp-serie-D2.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [53]:
fac = FAC(len(d2zpre), acf(d2zpre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos:
r12: -0.6150182129084292
r57: -0.10922275236044279
r69: 0.16716305731426437
r73: -0.11365881167509073
r81: -0.12386243150076603
r82: -0.11601226324982332
In [54]:
facp = FACP(len(d2zpre), pacf(d2zpre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos:
rho 12: -0.6294838987246432
rho 17: -0.09052048792417981
rho 24: -0.5103011546739746
rho 29: -0.093217025773514
rho 36: -0.3551238626796082
rho 48: -0.2751156800940471
rho 57: -0.108916168798713
rho 60: -0.2379289221214959
rho 61: 0.20095345338550902
rho 65: -0.0834129475389058
rho 72: -0.29986531752257994
rho 73: 0.302496719075776
rho 74: -0.15214692323832477
rho 75: 0.12158013142041356
rho 77: -0.21501383907290988
rho 79: -0.08587576253968264
rho 80: 0.10575443412592936
rho 81: -0.11295475178675196
rho 83: 0.1412351900166441
rho 84: -0.47389037358187136
rho 85: 0.8052856123151159
rho 86: -2.7039482046544805
rho 87: -1.3736483987654273
rho 88: 1.0019992843716
rho 89: 235.89044017004235
rho 90: -0.9971937437125793
rho 91: 0.48021486078669995
rho 93: -0.20196306057247534
rho 94: 0.28532937400900626
rho 95: -0.23402172704695934
rho 101: 0.25176068174227584
rho 102: -0.5167154340590121
rho 103: 0.7690418805288879
rho 104: -1.3401560972346027
rho 105: -2.166656395549073
rho 106: 1.1988394283045958
rho 107: 2.0506797927028866
rho 108: -1.2599924429780232
rho 109: -2.603184507263249
rho 110: 0.5135974713732451
rho 111: 0.2590435314099863
rho 113: -0.3628723780659208
rho 114: 0.4207095680305757
rho 115: -0.3885886584070735
rho 116: 0.2697703099821219
rho 119: -0.12676720968125227
rho 120: 0.16418782499731185

TERCERA DIFERENCIA ESTACIONAL¶

In [55]:
# Primera diferencia estacional con periodicidad 12
d3zpre = d2zpre.diff(12).dropna() 
In [56]:
fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(d3zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(d3zpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=3 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-04-fac-facp-serie-D3.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [57]:
fac = FAC(len(d3zpre), acf(d3zpre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos:
r12: -0.6989860750076662
r24: 0.19148304290921273
r57: -0.13077245780878802
r69: 0.18180783859402475
r81: -0.1446388376312966
In [58]:
facp = FACP(len(d3zpre), pacf(d3zpre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos:
rho 12: -0.7131085439918823
rho 24: -0.6165329897683706
rho 29: -0.08067623128823141
rho 36: -0.5171728973407944
rho 37: 0.1982055108755419
rho 40: 0.10015883076893528
rho 41: -0.09228473683010521
rho 47: 0.09394099253473005
rho 48: -0.4929991451837127
rho 49: 0.27281984306019774
rho 50: -0.09255185942814254
rho 51: -0.10940745619258534
rho 52: 0.34514660641128075
rho 53: -0.3969949656195508
rho 54: 0.33091157283957406
rho 55: -0.24725105044265083
rho 57: -0.1678397209757971
rho 58: -0.41066728720157253
rho 59: 1.4752062518993367
rho 60: 3.447013197287627
rho 61: -0.965945384029374
rho 62: 6.762979431106962
rho 63: 1.0827131685786548
rho 64: -0.8201163848967379
rho 65: 1.5200084948857941
rho 66: 1.3863162790341865
rho 67: -2.332119329595056
rho 68: -1.2962242829344859
rho 69: 0.5251991629964926
rho 70: 0.12535213223464997
rho 71: -0.43609447528833045
rho 72: 0.35806234464280917
rho 73: -0.1138999259876181
rho 74: -0.10235078255767406
rho 75: 0.31494875075818757
rho 76: -0.20483724449372892
rho 77: 0.0921074399984271
rho 78: -0.08928014769113522
rho 79: 0.0854845471954933
rho 80: -0.28927169831303046
rho 81: 0.20932339081585724
rho 82: -0.11854528891185322
rho 84: 0.19492618448545254
rho 87: 0.21562177268699964
rho 88: -0.19739536579003492
rho 89: 0.1375019109969488
rho 90: -0.21573215769368453
rho 91: 0.13649794688884304
rho 92: -0.4130856229291036
rho 93: 0.5769402061332968
rho 94: -1.1748169850119188
rho 95: -6.0784018693922
rho 96: 0.9057557611685417
rho 97: -0.5526621327578615
rho 98: 0.3652101944200118
rho 99: -0.14934082325742695
rho 100: -0.20792999327433961
rho 101: 0.32248566533747575
rho 102: -0.43340585213601085
rho 103: 0.35823263199853145
rho 104: -0.16127048470991623
rho 105: -0.09034092543015515
rho 106: 0.8633154494572077
rho 107: -9.435548190300198
rho 108: -1.0889219702203907
rho 109: 1.4801001973243877
rho 110: 2.1038789559878754
rho 111: -0.7905244792702399
rho 112: 0.6723557185188808
rho 113: -0.9378271953920853
rho 114: 5.770582312925257
rho 115: 0.947663800814546
rho 116: 3.6394621713184363
rho 117: -0.8146958649632807
rho 118: 0.18972318689108425
rho 119: 0.4988023348063774
rho 120: -0.48455031039030716

COMPARACIÓN DE VARIANZAS CON DISTINTAS DIFERENCIAS ESTACIONALES¶

In [59]:
# Desviación estándar de la serie
stds = [np.std(zpre), np.std(dzpre), np.std(d2zpre), np.std(d3zpre)]

# Plotear la desviación estándar
plt.figure()
plt.bar(['D=0', 'D=1', 'D=2', 'D=3'], stds, color=[sns.color_palette()[12], sns.color_palette()[10], sns.color_palette()[8], sns.color_palette()[7]])
plt.ylabel('Desviación Estándar')
plt.title('Desviación Estándar de las Series Estandarizadas, Diferenciadas Precipitación Mensual en Puebla')
plt.savefig('imagenes/03-05-std-Ds.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

Nos quedamos con la primera diferencia estacional, porque tiene una mejor forma su FAC y FACP, además de tener menor varianza.

SERIE CON TRANSFORMACIÓN¶

Definición matemática de la transformación Yeo-Johnson, como se encuentra en la literatura estadística (propuesta por Ingram Olkin y I. Paul Yeo & R. J. Johnson en 2000).

Sea $x \in \mathbb{R}$ un valor (puede ser negativo, cero o positivo). La transformación $T(x; \lambda)$ se define como:

$$ T(x; \lambda) = \begin{cases} \frac{[(x + 1)^\lambda - 1]}{\lambda} & \text{si } x \geq 0, \, \lambda \ne 0 \\ \log(x + 1) & \text{si } x \geq 0, \, \lambda = 0 \\ -\frac{[(-x + 1)^{2 - \lambda} - 1]}{2 - \lambda} & \text{si } x < 0, \, \lambda \ne 2 \\ -\log(-x + 1) & \text{si } x < 0, \, \lambda = 2 \end{cases} $$

ESTIMADOR PUNTUAL¶

In [60]:
from sklearn.preprocessing import PowerTransformer
import numpy as np
In [61]:
X = pre.values.reshape(-1, 1)

# Ajuste original para obtener lambda estimado
pt = PowerTransformer(method='yeo-johnson', standardize=True)
ypre = pt.fit_transform(X)
lambda_est = pt.lambdas_[0]
In [62]:
print("Lambda estimado:", pt.lambdas_)  # Obtener el valor lambda estimado
Lambda estimado: [0.23849441]

INTERVALO DE CONFIANZA¶

In [63]:
# Bootstrap para el IC de lambda
n_boot = 1000
lambdas_boot = []

for _ in range(n_boot):
    sample = np.random.choice(X.flatten(), size=len(X), replace=True).reshape(-1, 1)
    pt_boot = PowerTransformer(method='yeo-johnson', standardize=False)
    pt_boot.fit(sample)
    lambdas_boot.append(pt_boot.lambdas_[0])

# Calcular percentiles para el intervalo de confianza al 95%
ci_lower = np.percentile(lambdas_boot, 2.5)
ci_upper = np.percentile(lambdas_boot, 97.5)

print(f"Estimación de λ: {lambda_est:.4f}")
print(f"IC 95% para λ: ({ci_lower:.4f}, {ci_upper:.4f})")
Estimación de λ: 0.2385
IC 95% para λ: (0.2097, 0.2693)
In [64]:
# Gráfico opcional del histograma de los λ bootstrap
plt.figure(figsize=(10, 6))
plt.fill_betweenx([0, 120], ci_lower, ci_upper, color=sns.color_palette()[4], alpha=0.15)
plt.hist(lambdas_boot, bins=30, alpha=0.7, edgecolor='white')
plt.axvline(ci_lower, label=f'2.5% = {ci_lower:.4f}', color=sns.color_palette()[4], linewidth=3)
plt.axvline(ci_upper, label=f'97.5% = {ci_upper:.4f}', color=sns.color_palette()[4], linewidth=3)
plt.axvline(lambda_est, label=f'λ original = {lambda_est:.4f}', color=sns.color_palette()[1], linewidth=3)
plt.title('Distribución Bootstrap de λ')
plt.xlabel('λ')
plt.ylabel('Frecuencia')
plt.legend()
plt.gca().xaxis.set_major_locator(plt.MultipleLocator(0.01))
plt.savefig('imagenes/02-lambda-yj-bootstrapping.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

SERIE TRANSFORMADA¶

In [65]:
plt.figure()
plt.plot(date[:-12], ypre)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Yeo-Johnson Estandarizado Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-precipitacion_yeo_johnson.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

ESTACIONARIEDAD¶

In [66]:
from statsmodels.tsa.stattools import adfuller  
In [67]:
adfuller(ypre)
Out[67]:
(-6.885194656270816,
 1.4004074124320326e-09,
 16,
 643,
 {'1%': -3.440560883168159,
  '5%': -2.8660454146233434,
  '10%': -2.569169329058723},
 1139.5328800342036)
Estadístico ADF         = -6.885  
Valor-p                 = 1.40e-09  
Número de rezagos       = 16  
Número de observaciones = 643  
Valores críticos:
    1%  -> -3.4406
    5%  -> -2.8660
    10% -> -2.5692
Log-likelihood          = 3742.248
  • El estadístico ADF es $-6.885$, menor que los valores críticos a los niveles del 1%.
  • El valor-p es $1.40 × 10^{-9}$ mucho menor al nivel de significancia $\alpha$ = 0.05.

Por lo tanto, se rechaza la hipótesis nula, es decir, se concluye que la serie es estacionaria.

AUTOCORRELACIONES¶

In [ ]:
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf, acf, pacf
import matplotlib.pyplot as plt
from fac_y_facps_significativas import *    # Este es el archivo que se usa para el curso
In [ ]:
fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
axs[0].set_xlabel('Retraso')
axs[0].set_ylabel('ACF')

plot_pacf(ypre, lags=48, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].set_xlabel('Retraso')
axs[1].set_ylabel('PACF')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Transformada de la Precipitación Mensual en Puebla', fontsize=20, weight='bold')
plt.savefig('imagenes/04-acf_pacf.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [68]:
fac = FAC(len(ypre), acf(ypre, nlags=48)[1:] )
Valores de autocorrelacion significativos:
r1: 0.6331507478065659
r2: 0.3325589848089995
r4: -0.3429619382624403
r5: -0.6128788441981868
r6: -0.7155687723661789
r7: -0.5667143985229659
r8: -0.28172958031763784
r10: 0.34639100385471927
r11: 0.6139734232526982
r12: 0.6833003454277807
r13: 0.5238197525884694
r14: 0.23526370655336934
r16: -0.37665249496039976
r17: -0.6133232329123165
r18: -0.6522050881633413
r19: -0.49864844118107393
r22: 0.3581357207867556
r23: 0.6051041828734448
r24: 0.632952748472827
r25: 0.4787907060076014
r28: -0.3830626619493817
r29: -0.5771240646474324
r30: -0.6219569734501011
r31: -0.44116167635231773
r34: 0.3905997714697265
r35: 0.5842211455257492
r36: 0.606315712870404
r37: 0.41306326727063364
r40: -0.39053448733579516
r41: -0.5677702589671872
r42: -0.5689341833781766
r43: -0.3981104753267791
r46: 0.3759283566696746
r47: 0.5704749093905213
r48: 0.5490665087242427
In [69]:
facp = FACP(len(ypre), pacf(ypre, nlags=48)[1:])
Valores de autocorrelacion parcial significativos:
rho 1: 0.6341115228411734
rho 2: -0.11461335623878138
rho 3: -0.3000454427793004
rho 4: -0.33135817233932896
rho 5: -0.38109982701434264
rho 6: -0.32479269111279224
rho 7: -0.08061163300684875
rho 11: 0.2000309079008965
rho 12: 0.1532581642408276
rho 17: -0.12292322538615731
rho 22: -0.08358889406877636
rho 23: 0.12630684698989092
rho 30: -0.1028445923709159
rho 47: 0.09892532485023141

DIFERENCIAS ESTACIONALES¶

Por la FAC y la FACP, la serie transformada también necesita diferencias estacionales.

PRIMERA DIFERENCIA ESTACIONAL¶

In [71]:
# Primera diferencia estacional con periodicidad 12
ypre_series = pd.Series(ypre.flatten(), index=pre.index)       # Convertimos a serie de tiempo
dypre = ypre_series.diff(12).dropna()                          # Diferencia estacional
In [72]:
# Visualizamos la serie
dypre
Out[72]:
195209
1954-02-01   -0.109427
1954-03-01   -1.152261
1954-04-01    1.611746
1954-05-01    1.819524
1954-06-01    1.900692
                ...   
2008-08-01    0.606731
2008-09-01    0.156192
2008-10-01    0.052883
2008-11-01    0.312361
2008-12-01    0.000000
Length: 648, dtype: float64
In [74]:
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3))  # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha
plt.plot(dypre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Primera Diferencia Estacional, Transformación Estandarizada, Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-01-primera-diferencia-estacional-precipitacion.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [76]:
adfuller(dypre)
Out[76]:
(-9.030402747697376,
 5.464639337632109e-15,
 16,
 631,
 {'1%': -3.440755866431696,
  '5%': -2.86613130039063,
  '10%': -2.569215089800357},
 1308.464428619442)
Estadístico ADF         = -10.871  
Valor-p                 = 1.37e-19  
Número de rezagos       = 11  
Número de observaciones = 636  
Valores críticos:
    1%  -> -3.4407
    5%  -> -2.8661
    10% -> -2.5692
Log-likelihood          = 9952.186
  • El estadístico ADF es $-10.871$, menor que los valores críticos a los niveles del 1%.
  • El valor-p es $1.37 × 10^{-19}$ mucho menor al nivel de significancia $\alpha$ = 0.05.

Por lo tanto, se rechaza la hipótesis nula, es decir, se concluye que la serie es estacionaria.

In [77]:
fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(dypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(dypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=1, Estandarizada de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-02-fac-facp-serie-D1.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [78]:
fac = FAC(len(dypre), acf(dypre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos:
r1: 0.1881206416574384
r2: 0.12601374062799692
r6: -0.10023577930193152
r12: -0.40245032162240263
r33: 0.09802046981427195
r45: -0.11985655705785987
In [79]:
facp = FACP(len(dypre), pacf(dypre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos:
rho 1: 0.18841139999075743
rho 2: 0.09425081653205467
rho 6: -0.0848567418000625
rho 12: -0.4319892618507997
rho 13: 0.08136260875048439
rho 17: -0.11385209887045368
rho 22: -0.08401445717662205
rho 24: -0.2809725408252081
rho 25: 0.10755279636113212
rho 30: -0.09673228306484968
rho 33: 0.1105360450418002
rho 36: -0.14602967562386657
rho 42: -0.11769525648448687
rho 46: -0.09378255882994088
rho 48: -0.16710968828081296
rho 50: 0.09151867173223789
rho 55: 0.09278816335947565
rho 60: -0.09971581325023068
rho 72: -0.13002571917786468
rho 78: -0.07964448430687596
rho 84: -0.07986957194068202
rho 92: -0.13966165907823416
rho 93: 0.11600081821072342
rho 94: -0.1221857229746919
rho 113: -0.08611614522897519
rho 118: -0.08199830520693245
rho 120: -0.11627799482364788

SEGUNDA DIFERENCIA ESTACIONAL¶

In [80]:
# 2 diferencia estacional con periodicidad 12
d2ypre = dypre.diff(12).dropna() 
In [81]:
fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(d2ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(d2ypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')


plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=2 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-03-fac-facp-serie-D2.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [ ]:
fac = FAC(len(d2ypre), acf(d2ypre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos:
r12: -0.615018212908429
r57: -0.10922275236044277
r69: 0.16716305731426437
r73: -0.1136588116750907
r81: -0.12386243150076602
r82: -0.11601226324982332
In [ ]:
facp = FACP(len(d2ypre), pacf(d2ypre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos:
rho 12: -0.6294838987246429
rho 17: -0.09052048792417974
rho 24: -0.510301154673974
rho 29: -0.09321702577351394
rho 36: -0.35512386267960755
rho 48: -0.2751156800940455
rho 57: -0.10891616879871288
rho 60: -0.23792892212149455
rho 61: 0.20095345338550874
rho 65: -0.08341294753890602
rho 72: -0.29986531752257806
rho 73: 0.30249671907577647
rho 74: -0.1521469232383268
rho 75: 0.12158013142041249
rho 77: -0.21501383907291427
rho 79: -0.08587576253968643
rho 80: 0.10575443412592932
rho 81: -0.11295475178675143
rho 83: 0.1412351900166407
rho 84: -0.4738903735818693
rho 85: 0.8052856123151381
rho 86: -2.7039482046551293
rho 87: -1.373648398765204
rho 88: 1.0019992843716292
rho 89: 235.89044016524818
rho 90: -0.9971937437125115
rho 91: 0.4802148607866808
rho 93: -0.20196306057246763
rho 94: 0.2853293740090054
rho 95: -0.23402172704695937
rho 101: 0.25176068174227473
rho 102: -0.5167154340590162
rho 103: 0.769041880528836
rho 104: -1.3401560972344024
rho 105: -2.166656395548571
rho 106: 1.198839428304655
rho 107: 2.0506797927018146
rho 108: -1.259992442977829
rho 109: -2.6031845072646265
rho 110: 0.5135974713735127
rho 111: 0.2590435314096445
rho 113: -0.36287237806607014
rho 114: 0.4207095680306412
rho 115: -0.3885886584071046
rho 116: 0.26977030998213086
rho 119: -0.12676720968126257
rho 120: 0.16418782499731296

TERCERA DIFERENCIA ESTACIONAL¶

In [83]:
# Primera diferencia estacional con periodicidad 12
d3ypre = d2ypre.diff(12).dropna() 
In [84]:
fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(d3ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(d3ypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=3 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-04-fac-facp-serie-D3.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [85]:
fac = FAC(len(d3ypre), acf(d3ypre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos:
r1: 0.14142525988694732
r11: -0.08618353748437811
r12: -0.7006965212651104
r13: -0.12299455598817151
r24: 0.21088909322401359
r33: 0.14799260300636627
r34: 0.124164682223292
r45: -0.13840321103959075
In [86]:
facp = FACP(len(d3ypre), pacf(d3ypre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos:
rho 1: 0.14165226672464704
rho 11: -0.08572854737384034
rho 12: -0.7145741931403762
rho 14: 0.10244937563798918
rho 21: -0.0974412193826169
rho 22: -0.11958849854177123
rho 24: -0.590999784962588
rho 25: 0.1155406109341616
rho 36: -0.42596863520035666
rho 37: 0.18734873402939595
rho 42: -0.1733829621796275
rho 46: -0.1132553928924406
rho 48: -0.39708277000220765
rho 49: 0.17397913808655488
rho 54: -0.3269257232330364
rho 55: 0.21559931119807146
rho 56: -0.19500988553521115
rho 58: -0.22348004365902227
rho 59: 0.25219826387781036
rho 60: -0.7375993663060407
rho 61: 2.2537091363236095
rho 62: 1.650176227081085
rho 63: -0.6968349029252933
rho 64: 0.3967431321320999
rho 65: -0.2123951680489383
rho 66: -0.20497281303308482
rho 68: 0.23621005007258603
rho 69: -0.8189484098245967
rho 70: 1.3011193720883465
rho 71: 1.4316135793533165
rho 72: -5.898535785106428
rho 73: -1.3199340565173563
rho 74: -0.6772780595533806
rho 75: 0.5245892189329673
rho 76: -0.2815596779482279
rho 77: -0.4521464433566569
rho 78: 0.30058897281376856
rho 79: -0.4164535343317919
rho 80: -0.11859502105379839
rho 81: -0.15366212304252277
rho 82: -0.2922622187005957
rho 83: -0.4928118125453261
rho 84: -0.14977995949488318
rho 85: -0.8927516420039372
rho 86: -5.928900128261838
rho 87: 1.2117900757453872
rho 88: 0.5145451003205691
rho 89: 0.24638896698250676
rho 90: 0.5283271371105461
rho 92: 0.3314944963217833
rho 93: 0.11506584057595354
rho 94: 0.21220949431027086
rho 95: -0.12891904223431166
rho 96: 0.5141635439947064
rho 97: -0.6188892690196736
rho 98: 1.7101159305305729
rho 99: 2.12107036942252
rho 100: -0.6062698236020745
rho 101: 0.46103150897488143
rho 102: -0.1569944090737856
rho 103: 0.2535125462534724
rho 104: -0.21007258043287413
rho 105: 0.3130566565661735
rho 107: 0.1210008208606222
rho 109: 0.1042542157962568
rho 110: -0.17221115272098989
rho 111: 0.2646586117130243
rho 112: -0.11554768465356624
rho 114: 0.18496589221568
rho 115: 0.10819324030267996
rho 116: -0.15857098929802355
rho 117: 0.3781749655329778
rho 118: -0.1243216624045143
rho 120: 0.20596573821321837

COMPARACIÓN DE VARIANZAS CON DISTINTAS DIFERENCIAS ESTACIONALES¶

In [88]:
# Desviación estándar de la serie
stds = [np.std(ypre), np.std(dypre), np.std(d2ypre), np.std(d3ypre)]

# Plotear la desviación estándar
plt.figure()
plt.bar(['D=0', 'D=1', 'D=2', 'D=3'], stds, color=[sns.color_palette()[12], sns.color_palette()[10], sns.color_palette()[8], sns.color_palette()[7]])
plt.ylabel('Desviación Estándar')
plt.title('Desviación Estándar de las Series Estandarizadas, Diferenciadas Precipitación Mensual en Puebla')
plt.savefig('imagenes/03-05-std-Ds.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

Nos quedamos con la primera diferencia estacional, porque tiene una mejor forma su FAC y FACP, además de tener menor varianza.

MODELADO¶

Lo que se muestra a continuación fue realizado para la serie sin transformar, solamente estandarizada.

En 02-----Encontrar-los-Mejores-Modelos se utilizá la libreria de pmdarima para intentar diferentes combinaciones de modelos y encontrar sus $\text{AIC}$.

Se obtuvieron los siguientes resultados:

(1)  ARIMA(1,0,0)(0,1,1)[12], AIC=1430.397
(2)  ARIMA(0,0,1)(0,1,1)[12], AIC=1431.266
(3)  ARIMA(1,0,1)(0,1,1)[12], AIC=1432.133
(4)  ARIMA(2,0,0)(0,1,1)[12], AIC=1432.157
(5)  ARIMA(1,0,0)(0,1,2)[12], AIC=1432.368
(6)  ARIMA(1,0,0)(1,1,1)[12], AIC=1432.373
(7)  ARIMA(0,0,2)(0,1,1)[12], AIC=1432.664
(8)  ARIMA(0,0,1)(0,1,2)[12], AIC=1433.240
(9)  ARIMA(0,0,1)(1,1,1)[12], AIC=1433.245
(10) ARIMA(2,0,1)(0,1,1)[12], AIC=1434.127
(11) ARIMA(0,0,0)(0,1,1)[12], AIC=1445.334
(12) ARIMA(1,0,0)(1,1,0)[12], AIC=1523.141
(13) ARIMA(0,0,1)(1,1,0)[12], AIC=1523.219
(14) ARIMA(0,0,1)(0,1,0)[12], AIC=1652.849
(15) ARIMA(1,0,0)(0,1,0)[12], AIC=1652.984
(16) ARIMA(0,0,0)(0,1,0)[12], AIC=1655.671

Ningun modelo cumplió los supuestos. 02-01-Encontrar-Modelos-que-Cumplan

Por la FAC y FACP, se creyó que podía ser un $\text{ARIMA}(1, 0, 1) \times (0, 1, 1, 12)$, pero este modelo, con la prueba de Ljung-Box, tenía dependencia entre los residuos. Ver 03-----Encontrar-Modelos-que-Cumplan, 1 MODELO.

Se propuso aumentar el número de parámetros, se pusieron parámetros "intuitivos", es decir, la precipitación depende de los últimos 5 meses, de los últimos 4 años, y hay residuos en los últimos dos meses, llegando al primer modelo que cumplía la independecia en los RESIDUOS siendo el $\text{ARIMA}(5, 0, 2) \times (4, 1, 0, 12)$. Sin embargo, no era estacionario ni invertible. Ver 03-----Encontrar-Modelos-que-Cumplan, 2 MODELO.

Se fue variando el número de parámetros de los modelos, notando curiosamente, que para esta serie, que los modelos que eran estacionarios eran los que tenían su parte del polinomio autorregresivo completa y no tenían parte autorregresiva estacional (y viceversa), lo mismo para la parte media móvil y la invertibilidad. Tomando en cuenta esta observación, se hicieron las notebooks de 03-01, ..., 03-04, donde se varían, por ejemplo, los $p$, manteniendo $P=0$, lo mismo para la parte media móvil. Lo malo, es que para los modelos tuvieran independencia en los residuos, los polinomios tenían que ser bastante grandes, lo que podía hacer que no se cumpliera el principio de parsimonia.

En las notebooks 03-01, ..., 03-04 se proponen multiples modelos, que, dada la Observación (2), que dice que para esta serie los modelos que cumplen esa forma son admisibles, se siguió el siguiente procedimiento:

  • Se fijan los valores de $\text{ARIMA}(p, 0, q) \times (0, 1, 0, 12)$ donde $p = 0, 1, 2, \dots$ y $q = 0, 1, 2, \dots$.
  • Se fijan los valores de $\text{ARIMA}(0, 0, q) \times (P, 1, 0, 12)$ donde $P = 0, 1, 2, \dots$ y $q = 0, 1, 2, \dots$.
  • Se fijan los valores de $\text{ARIMA}(p, 0, 0) \times (0, 1, Q, 12)$ donde $p = 0, 1, 2, \dots$ y $Q = 0, 1, 2, \dots$.
  • Se fijan los valores de $\text{ARIMA}(0, 0, 0) \times (P, 1, Q, 12)$ donde $P = 0, 1, 2, \dots$ y $Q = 0, 1, 2, \dots$.

RESULTADOS de 03-01-Encontrar-Modelo-pq:

(1)   ARIMA(0,0,15)(0,1,0)[12]    AIC=1370.81  
(2)   ARIMA(0,0,12)(0,1,0)[12]    AIC=1370.92  
(3)   ARIMA(3,0,12)(0,1,0)[12]    AIC=1372.25  
(4)   ARIMA(0,0,16)(0,1,0)[12]    AIC=1372.52  
(5)   ARIMA(1,0,15)(0,1,0)[12]    AIC=1372.60  
(6)   ARIMA(1,0,12)(0,1,0)[12]    AIC=1372.69  
(7)   ARIMA(0,0,13)(0,1,0)[12]    AIC=1372.78  
(8)   ARIMA(3,0,13)(0,1,0)[12]    AIC=1372.87  
(9)   ARIMA(4,0,12)(0,1,0)[12]    AIC=1373.41  
(10)  ARIMA(1,0,16)(0,1,0)[12]    AIC=1373.59  
(11)  ARIMA(1,0,13)(0,1,0)[12]    AIC=1373.76  
(12)  ARIMA(3,0,14)(0,1,0)[12]    AIC=1374.01  
(13)  ARIMA(4,0,13)(0,1,0)[12]    AIC=1374.12  
(14)  ARIMA(5,0,12)(0,1,0)[12]    AIC=1374.34  
(15)  ARIMA(0,0,14)(0,1,0)[12]    AIC=1374.44  
(16)  ARIMA(2,0,12)(0,1,0)[12]    AIC=1374.49  
(17)  ARIMA(2,0,12)(0,1,0)[12]    AIC=1374.49  
(18)  ARIMA(0,0,17)(0,1,0)[12]    AIC=1374.49  
(19)  ARIMA(2,0,15)(0,1,0)[12]    AIC=1374.62  
(20)  ARIMA(2,0,15)(0,1,0)[12]    AIC=1374.62  
(21)  ARIMA(1,0,17)(0,1,0)[12]    AIC=1374.97  
(22)  ARIMA(4,0,14)(0,1,0)[12]    AIC=1375.08  
(23)  ARIMA(2,0,13)(0,1,0)[12]    AIC=1375.16  
(24)  ARIMA(2,0,13)(0,1,0)[12]    AIC=1375.16  
(25)  ARIMA(1,0,14)(0,1,0)[12]    AIC=1375.26  
(26)  ARIMA(5,0,13)(0,1,0)[12]    AIC=1375.43  
(27)  ARIMA(6,0,12)(0,1,0)[12]    AIC=1375.94  
(28)  ARIMA(6,0,12)(0,1,0)[12]    AIC=1375.94  
(29)  ARIMA(2,0,16)(0,1,0)[12]    AIC=1375.97  
(30)  ARIMA(2,0,16)(0,1,0)[12]    AIC=1375.97  
(31)  ARIMA(3,0,15)(0,1,0)[12]    AIC=1376.37  
(32)  ARIMA(7,0,12)(0,1,0)[12]    AIC=1376.46  
(33)  ARIMA(2,0,14)(0,1,0)[12]    AIC=1376.63  
(34)  ARIMA(2,0,14)(0,1,0)[12]    AIC=1376.63  
(35)  ARIMA(5,0,14)(0,1,0)[12]    AIC=1377.04  
(36)  ARIMA(6,0,13)(0,1,0)[12]    AIC=1377.52  
(37)  ARIMA(6,0,13)(0,1,0)[12]    AIC=1377.52  
(38)  ARIMA(4,0,15)(0,1,0)[12]    AIC=1377.61  
(39)  ARIMA(2,0,17)(0,1,0)[12]    AIC=1377.67  
(40)  ARIMA(2,0,17)(0,1,0)[12]    AIC=1377.67  
(41)  ARIMA(3,0,16)(0,1,0)[12]    AIC=1377.71  
(42)  ARIMA(7,0,13)(0,1,0)[12]    AIC=1378.15  
(43)  ARIMA(8,0,12)(0,1,0)[12]    AIC=1378.47  
(44)  ARIMA(8,0,12)(0,1,0)[12]    AIC=1378.47  
(45)  ARIMA(9,0,12)(0,1,0)[12]    AIC=1378.79  
(46)  ARIMA(6,0,14)(0,1,0)[12]    AIC=1378.97  
(47)  ARIMA(6,0,14)(0,1,0)[12]    AIC=1378.97  
(48)  ARIMA(4,0,16)(0,1,0)[12]    AIC=1379.46  
(49)  ARIMA(3,0,17)(0,1,0)[12]    AIC=1379.51  
(50)  ARIMA(5,0,15)(0,1,0)[12]    AIC=1379.68  
(51)  ARIMA(8,0,13)(0,1,0)[12]    AIC=1380.15  
(52)  ARIMA(8,0,13)(0,1,0)[12]    AIC=1380.15  
(53)  ARIMA(7,0,14)(0,1,0)[12]    AIC=1380.53  
(54)  ARIMA(10,0,12)(0,1,0)[12]   AIC=1380.62  
(55)  ARIMA(5,0,16)(0,1,0)[12]    AIC=1380.90  
(56)  ARIMA(4,0,17)(0,1,0)[12]    AIC=1381.04  
(57)  ARIMA(6,0,15)(0,1,0)[12]    AIC=1381.49  
(58)  ARIMA(6,0,15)(0,1,0)[12]    AIC=1381.49  
(59)  ARIMA(10,0,13)(0,1,0)[12]   AIC=1381.66  
(60)  ARIMA(9,0,13)(0,1,0)[12]    AIC=1381.75  
(61)  ARIMA(8,0,14)(0,1,0)[12]    AIC=1382.31  
(62)  ARIMA(8,0,14)(0,1,0)[12]    AIC=1382.31  
(63)  ARIMA(5,0,17)(0,1,0)[12]    AIC=1382.35  
(64)  ARIMA(11,0,12)(0,1,0)[12]   AIC=1382.48  
(65)  ARIMA(9,0,14)(0,1,0)[12]    AIC=1382.58  
(66)  ARIMA(6,0,16)(0,1,0)[12]    AIC=1382.64  
(67)  ARIMA(6,0,16)(0,1,0)[12]    AIC=1382.64  
(68)  ARIMA(7,0,15)(0,1,0)[12]    AIC=1382.73  
(69)  ARIMA(10,0,14)(0,1,0)[12]   AIC=1383.53  
(70)  ARIMA(11,0,13)(0,1,0)[12]   AIC=1383.82  
(71)  ARIMA(6,0,17)(0,1,0)[12]    AIC=1384.31  
(72)  ARIMA(6,0,17)(0,1,0)[12]    AIC=1384.31  
(73)  ARIMA(7,0,16)(0,1,0)[12]    AIC=1384.33  
(74)  ARIMA(8,0,15)(0,1,0)[12]    AIC=1384.52  
(75)  ARIMA(8,0,15)(0,1,0)[12]    AIC=1384.52  
(76)  ARIMA(9,0,15)(0,1,0)[12]    AIC=1384.53  
(77)  ARIMA(13,0,12)(0,1,0)[12]   AIC=1384.56  
(78)  ARIMA(11,0,14)(0,1,0)[12]   AIC=1384.73  
(79)  ARIMA(13,0,13)(0,1,0)[12]   AIC=1385.12  
(80)  ARIMA(10,0,15)(0,1,0)[12]   AIC=1385.45  
(81)  ARIMA(12,0,13)(0,1,0)[12]   AIC=1385.52  
(82)  ARIMA(8,0,16)(0,1,0)[12]    AIC=1385.88  
(83)  ARIMA(8,0,16)(0,1,0)[12]    AIC=1385.88  
(84)  ARIMA(12,0,12)(0,1,0)[12]   AIC=1386.15  
(85)  ARIMA(7,0,17)(0,1,0)[12]    AIC=1386.37  
(86)  ARIMA(11,0,15)(0,1,0)[12]   AIC=1386.42  
(87)  ARIMA(9,0,16)(0,1,0)[12]    AIC=1386.43  
(88)  ARIMA(13,0,14)(0,1,0)[12]   AIC=1386.89  
(89)  ARIMA(12,0,14)(0,1,0)[12]   AIC=1386.90  
(90)  ARIMA(10,0,16)(0,1,0)[12]   AIC=1387.11  
(91)  ARIMA(14,0,13)(0,1,0)[12]   AIC=1387.36  
(92)  ARIMA(14,0,12)(0,1,0)[12]   AIC=1387.85  
(93)  ARIMA(11,0,16)(0,1,0)[12]   AIC=1387.97  
(94)  ARIMA(8,0,17)(0,1,0)[12]    AIC=1388.02  
(95)  ARIMA(8,0,17)(0,1,0)[12]    AIC=1388.02  
(96)  ARIMA(15,0,12)(0,1,0)[12]   AIC=1388.42  
(97)  ARIMA(13,0,15)(0,1,0)[12]   AIC=1388.76  
(98)  ARIMA(12,0,15)(0,1,0)[12]   AIC=1388.84  
(99)  ARIMA(14,0,14)(0,1,0)[12]   AIC=1388.89  
(100) ARIMA(9,0,17)(0,1,0)[12]    AIC=1388.95  
(101) ARIMA(10,0,17)(0,1,0)[12]   AIC=1389.04  
(102) ARIMA(15,0,13)(0,1,0)[12]   AIC=1389.30  
(103) ARIMA(17,0,12)(0,1,0)[12]   AIC=1389.57  
(104) ARIMA(11,0,17)(0,1,0)[12]   AIC=1389.72  
(105) ARIMA(16,0,12)(0,1,0)[12]   AIC=1389.87  
(106) ARIMA(12,0,16)(0,1,0)[12]   AIC=1390.33  
(107) ARIMA(15,0,14)(0,1,0)[12]   AIC=1390.73  
(108) ARIMA(14,0,15)(0,1,0)[12]   AIC=1390.78  
(109) ARIMA(17,0,13)(0,1,0)[12]   AIC=1390.81  
(110) ARIMA(16,0,13)(0,1,0)[12]   AIC=1390.83  
(111) ARIMA(13,0,16)(0,1,0)[12]   AIC=1391.95  
(112) ARIMA(12,0,17)(0,1,0)[12]   AIC=1392.26  
(113) ARIMA(13,0,17)(0,1,0)[12]   AIC=1392.28  
(114) ARIMA(17,0,14)(0,1,0)[12]   AIC=1392.49  
(115) ARIMA(16,0,14)(0,1,0)[12]   AIC=1392.69  
(116) ARIMA(15,0,15)(0,1,0)[12]   AIC=1392.85  
(117) ARIMA(14,0,16)(0,1,0)[12]   AIC=1392.97  
(118) ARIMA(16,0,15)(0,1,0)[12]   AIC=1394.72  
(119) ARIMA(17,0,15)(0,1,0)[12]   AIC=1394.84  
(120) ARIMA(15,0,16)(0,1,0)[12]   AIC=1395.10  
(121) ARIMA(14,0,17)(0,1,0)[12]   AIC=1395.38  
(122) ARIMA(15,0,17)(0,1,0)[12]   AIC=1396.22  
(123) ARIMA(17,0,16)(0,1,0)[12]   AIC=1396.75  
(124) ARIMA(16,0,16)(0,1,0)[12]   AIC=1397.15  
(125) ARIMA(16,0,17)(0,1,0)[12]   AIC=1398.28  
(126) ARIMA(17,0,17)(0,1,0)[12]   AIC=1399.30  
(127) ARIMA(9,0,11)(0,1,0)[12]    AIC=1405.22  
(128) ARIMA(11,0,10)(0,1,0)[12]   AIC=1405.63  
(129) ARIMA(11,0,11)(0,1,0)[12]   AIC=1406.36  
(130) ARIMA(10,0,11)(0,1,0)[12]   AIC=1411.83  
(131) ARIMA(14,0,11)(0,1,0)[12]   AIC=1413.61  
(132) ARIMA(13,0,11)(0,1,0)[12]   AIC=1415.15  
(133) ARIMA(8,0,11)(0,1,0)[12]    AIC=1415.83  
(134) ARIMA(8,0,11)(0,1,0)[12]    AIC=1415.83  
(135) ARIMA(15,0,11)(0,1,0)[12]   AIC=1418.52  
(136) ARIMA(16,0,11)(0,1,0)[12]   AIC=1418.96  
(137) ARIMA(12,0,11)(0,1,0)[12]   AIC=1419.07  
(138) ARIMA(7,0,11)(0,1,0)[12]    AIC=1420.16  
(139) ARIMA(15,0,10)(0,1,0)[12]   AIC=1427.03  
(140) ARIMA(16,0,10)(0,1,0)[12]   AIC=1428.10  
(141) ARIMA(16,0,10)(0,1,0)[12]   AIC=1428.10  
(142) ARIMA(14,0,10)(0,1,0)[12]   AIC=1429.08  
(143) ARIMA(4,0,11)(0,1,0)[12]    AIC=1430.74  
(144) ARIMA(3,0,11)(0,1,0)[12]    AIC=1435.07  
(145) ARIMA(17,0,10)(0,1,0)[12]   AIC=1436.99  
(146) ARIMA(12,0,10)(0,1,0)[12]   AIC=1437.18  
(147) ARIMA(5,0,11)(0,1,0)[12]    AIC=1443.88  
(148) ARIMA(13,0,10)(0,1,0)[12]   AIC=1444.13  
(149) ARIMA(6,0,11)(0,1,0)[12]    AIC=1449.08  
(150) ARIMA(6,0,11)(0,1,0)[12]    AIC=1449.08  
(151) ARIMA(16,0,9)(0,1,0)[12]    AIC=1450.39  
(152) ARIMA(17,0,9)(0,1,0)[12]    AIC=1450.46  
(153) ARIMA(2,0,11)(0,1,0)[12]    AIC=1451.97  
(154) ARIMA(2,0,11)(0,1,0)[12]    AIC=1451.97  
(155) ARIMA(16,0,8)(0,1,0)[12]    AIC=1460.41  
(156) ARIMA(17,0,8)(0,1,0)[12]    AIC=1464.61  
(157) ARIMA(16,0,7)(0,1,0)[12]    AIC=1469.14  
(158) ARIMA(17,0,7)(0,1,0)[12]    AIC=1469.22  
(159) ARIMA(17,0,6)(0,1,0)[12]    AIC=1471.02  
(160) ARIMA(9,0,10)(0,1,0)[12]    AIC=1471.14  
(161) ARIMA(8,0,10)(0,1,0)[12]    AIC=1478.49  
(162) ARIMA(8,0,10)(0,1,0)[12]    AIC=1478.49  
(163) ARIMA(16,0,5)(0,1,0)[12]    AIC=1479.55  
(164) ARIMA(16,0,6)(0,1,0)[12]    AIC=1480.98  
(165) ARIMA(10,0,10)(0,1,0)[12]   AIC=1495.41  
(166) ARIMA(17,0,5)(0,1,0)[12]    AIC=1496.61  
(167) ARIMA(6,0,10)(0,1,0)[12]    AIC=1497.53  
(168) ARIMA(6,0,10)(0,1,0)[12]    AIC=1497.53  
(169) ARIMA(1,0,11)(0,1,0)[12]    AIC=1503.08  
(170) ARIMA(5,0,10)(0,1,0)[12]    AIC=1526.92  
(171) ARIMA(10,0,6)(0,1,0)[12]    AIC=1541.36  
(172) ARIMA(4,0,10)(0,1,0)[12]    AIC=1542.28  
(173) ARIMA(7,0,10)(0,1,0)[12]    AIC=1546.44  
(174) ARIMA(3,0,10)(0,1,0)[12]    AIC=1552.06  
(175) ARIMA(2,0,10)(0,1,0)[12]    AIC=1552.72  
(176) ARIMA(2,0,10)(0,1,0)[12]    AIC=1552.72  
(177) ARIMA(2,0,10)(0,1,0)[12]    AIC=1552.72  
(178) ARIMA(2,0,9)(0,1,0)[12]     AIC=1555.97  
(179) ARIMA(2,0,8)(0,1,0)[12]     AIC=1559.09  
(180) ARIMA(1,0,8)(0,1,0)[12]     AIC=1593.01  
(181) ARIMA(0,0,11)(0,1,0)[12]    AIC=1599.58  
(182) ARIMA(2,0,5)(0,1,0)[12]     AIC=1601.11  
(183) ARIMA(2,0,6)(0,1,0)[12]     AIC=1621.74  
(184) ARIMA(0,0,9)(0,1,0)[12]     AIC=1621.79  
(185) ARIMA(0,0,10)(0,1,0)[12]    AIC=1622.73  
(186) ARIMA(0,0,10)(0,1,0)[12]    AIC=1622.73  
(187) ARIMA(1,0,9)(0,1,0)[12]     AIC=1623.45  
(188) ARIMA(1,0,10)(0,1,0)[12]    AIC=1624.20  
(189) ARIMA(1,0,10)(0,1,0)[12]    AIC=1624.20  
(190) ARIMA(1,0,5)(0,1,0)[12]     AIC=1627.56  
(191) ARIMA(1,0,6)(0,1,0)[12]     AIC=1630.49  
(192) ARIMA(0,0,7)(0,1,0)[12]     AIC=1630.59  
(193) ARIMA(0,0,8)(0,1,0)[12]     AIC=1630.74  
(194) ARIMA(1,0,7)(0,1,0)[12]     AIC=1631.92  
(195) ARIMA(2,0,7)(0,1,0)[12]     AIC=1633.08  
(196) ARIMA(0,0,6)(0,1,0)[12]     AIC=1635.97  
(197) ARIMA(0,0,5)(0,1,0)[12]     AIC=1651.91

RESULTADOS de 03-02-Encontrar-Modelo-Pq:

(1)   ARIMA(0,0,15)(0,1,0)[12]    AIC=1370.81  
(2)   ARIMA(0,0,12)(0,1,0)[12]    AIC=1370.92  
(3)   ARIMA(0,0,12)(1,1,0)[12]    AIC=1372.79  
(4)   ARIMA(0,0,15)(1,1,0)[12]    AIC=1373.17  
(5)   ARIMA(0,0,12)(2,1,0)[12]    AIC=1373.37  
(6)   ARIMA(0,0,13)(1,1,0)[12]    AIC=1373.86  
(7)   ARIMA(0,0,12)(4,1,0)[12]    AIC=1375.93  
(8)   ARIMA(0,0,14)(2,1,0)[12]    AIC=1377.49  
(9)   ARIMA(0,0,11)(3,1,0)[12]    AIC=1427.35  
(10)  ARIMA(0,0,10)(3,1,0)[12]    AIC=1430.51  
(11)  ARIMA(0,0,11)(2,1,0)[12]    AIC=1443.12  
(12)  ARIMA(0,0,10)(2,1,0)[12]    AIC=1446.97  
(13)  ARIMA(0,0,11)(1,1,0)[12]    AIC=1511.21  
(14)  ARIMA(0,0,6)(1,1,0)[12]     AIC=1514.47  
(15)  ARIMA(0,0,8)(1,1,0)[12]     AIC=1516.15  
(16)  ARIMA(0,0,10)(1,1,0)[12]    AIC=1516.98

RESULTADOS de 03-03-Encontrar-Modelo-pQ:

(1)   ARIMA(12,0,0)(0,1,4)[12]    AIC=1367.04  
(2)   ARIMA(13,0,0)(0,1,4)[12]    AIC=1367.37  
(3)   ARIMA(12,0,0)(0,1,3)[12]    AIC=1368.24  
(4)   ARIMA(12,0,0)(0,1,5)[12]    AIC=1369.25  
(5)   ARIMA(15,0,0)(0,1,4)[12]    AIC=1371.49  
(6)   ARIMA(12,0,0)(0,1,1)[12]    AIC=1382.66  
(7)   ARIMA(15,0,0)(0,1,1)[12]    AIC=1387.68  
(8)   ARIMA(6,0,0)(0,1,1)[12]     AIC=1391.61  
(9)   ARIMA(11,0,0)(0,1,4)[12]    AIC=1392.63  
(10)  ARIMA(10,0,0)(0,1,4)[12]    AIC=1398.90  
(11)  ARIMA(3,0,0)(0,1,1)[12]     AIC=1431.92  
(12)  ARIMA(0,0,0)(0,1,1)[12]     AIC=1443.59  
(13)  ARIMA(12,0,0)(0,1,0)[12]    AIC=1522.91  
(14)  ARIMA(6,0,0)(0,1,0)[12]     AIC=1650.81

RESULTADOS de 03-04-Encontrar-Modelo-PQ:

(1)   ARIMA(0,0,0)(0,1,3)[12]     AIC=1443.50  
(2)   ARIMA(0,0,0)(0,1,1)[12]     AIC=1443.59  
(3)   ARIMA(0,0,0)(2,1,1)[12]     AIC=1444.11  
(4)   ARIMA(0,0,0)(1,1,2)[12]     AIC=1445.15  
(5)   ARIMA(0,0,0)(1,1,3)[12]     AIC=1445.28  
(6)   ARIMA(0,0,0)(0,1,4)[12]     AIC=1445.49  
(7)   ARIMA(0,0,0)(0,1,2)[12]     AIC=1445.59  
(8)   ARIMA(0,0,0)(1,1,1)[12]     AIC=1445.59  
(9)   ARIMA(0,0,0)(4,1,3)[12]     AIC=1445.70  
(10)  ARIMA(0,0,0)(4,1,4)[12]     AIC=1445.74  
(11)  ARIMA(0,0,0)(2,1,2)[12]     AIC=1445.92  
(12)  ARIMA(0,0,0)(4,1,2)[12]     AIC=1445.95  
(13)  ARIMA(0,0,0)(3,1,1)[12]     AIC=1446.07  
(14)  ARIMA(0,0,0)(2,1,3)[12]     AIC=1446.38  
(15)  ARIMA(0,0,0)(1,1,4)[12]     AIC=1446.57  
(16)  ARIMA(0,0,0)(4,1,1)[12]     AIC=1446.70  
(17)  ARIMA(0,0,0)(3,1,2)[12]     AIC=1446.99  
(18)  ARIMA(0,0,0)(3,1,3)[12]     AIC=1447.78  
(19)  ARIMA(0,0,0)(2,1,4)[12]     AIC=1447.80  
(20)  ARIMA(0,0,0)(4,1,0)[12]     AIC=1448.24  
(21)  ARIMA(0,0,0)(3,1,4)[12]     AIC=1449.21  
(22)  ARIMA(0,0,0)(3,1,0)[12]     AIC=1454.13  
(23)  ARIMA(0,0,0)(2,1,0)[12]     AIC=1464.14  
(24)  ARIMA(0,0,0)(1,1,0)[12]     AIC=1526.64  
(25)  ARIMA(0,0,0)(0,1,0)[12]     AIC=1653.67

Estos son todos los modelos propuestos, ordenados por su AIC:

(1) ARIMA(12,0,0)(0,1,4)[12], AIC=1367.04
(2) ARIMA(13,0,0)(0,1,4)[12], AIC=1367.37
(3) ARIMA(12,0,0)(0,1,3)[12], AIC=1368.24
(4) ARIMA(12,0,0)(0,1,5)[12], AIC=1369.25
(5) ARIMA(0,0,15)(0,1,0)[12], AIC=1370.81
(6) ARIMA(0,0,15)(0,1,0)[12], AIC=1370.81
(7) ARIMA(0,0,12)(0,1,0)[12], AIC=1370.92
(8) ARIMA(0,0,12)(0,1,0)[12], AIC=1370.92
(9) ARIMA(15,0,0)(0,1,4)[12], AIC=1371.49
(10) ARIMA(3,0,12)(0,1,0)[12], AIC=1372.25
(11) ARIMA(0,0,16)(0,1,0)[12], AIC=1372.52
(12) ARIMA(1,0,15)(0,1,0)[12], AIC=1372.60
(13) ARIMA(1,0,12)(0,1,0)[12], AIC=1372.69
(14) ARIMA(0,0,13)(0,1,0)[12], AIC=1372.78
(15) ARIMA(0,0,12)(1,1,0)[12], AIC=1372.79
(16) ARIMA(3,0,13)(0,1,0)[12], AIC=1372.87
(17) ARIMA(0,0,15)(1,1,0)[12], AIC=1373.17
(18) ARIMA(0,0,12)(2,1,0)[12], AIC=1373.37
(19) ARIMA(4,0,12)(0,1,0)[12], AIC=1373.41
(20) ARIMA(1,0,16)(0,1,0)[12], AIC=1373.59
(21) ARIMA(1,0,13)(0,1,0)[12], AIC=1373.76
(22) ARIMA(0,0,13)(1,1,0)[12], AIC=1373.86
(23) ARIMA(3,0,14)(0,1,0)[12], AIC=1374.01
(24) ARIMA(4,0,13)(0,1,0)[12], AIC=1374.12
(25) ARIMA(5,0,12)(0,1,0)[12], AIC=1374.34
(26) ARIMA(0,0,14)(0,1,0)[12], AIC=1374.44
(27) ARIMA(2,0,12)(0,1,0)[12], AIC=1374.49
(28) ARIMA(2,0,12)(0,1,0)[12], AIC=1374.49
(29) ARIMA(0,0,17)(0,1,0)[12], AIC=1374.49
(30) ARIMA(2,0,15)(0,1,0)[12], AIC=1374.62
(31) ARIMA(2,0,15)(0,1,0)[12], AIC=1374.62
(32) ARIMA(1,0,17)(0,1,0)[12], AIC=1374.97
(33) ARIMA(4,0,14)(0,1,0)[12], AIC=1375.08
(34) ARIMA(2,0,13)(0,1,0)[12], AIC=1375.16
(35) ARIMA(2,0,13)(0,1,0)[12], AIC=1375.16
(36) ARIMA(1,0,14)(0,1,0)[12], AIC=1375.26
(37) ARIMA(5,0,13)(0,1,0)[12], AIC=1375.43
(38) ARIMA(0,0,12)(4,1,0)[12], AIC=1375.93
(39) ARIMA(6,0,12)(0,1,0)[12], AIC=1375.94
(40) ARIMA(6,0,12)(0,1,0)[12], AIC=1375.94
(41) ARIMA(2,0,16)(0,1,0)[12], AIC=1375.97
(42) ARIMA(2,0,16)(0,1,0)[12], AIC=1375.97
(43) ARIMA(3,0,15)(0,1,0)[12], AIC=1376.37
(44) ARIMA(7,0,12)(0,1,0)[12], AIC=1376.46
(45) ARIMA(2,0,14)(0,1,0)[12], AIC=1376.63
(46) ARIMA(2,0,14)(0,1,0)[12], AIC=1376.63
(47) ARIMA(5,0,14)(0,1,0)[12], AIC=1377.04
(48) ARIMA(0,0,14)(2,1,0)[12], AIC=1377.49
(49) ARIMA(6,0,13)(0,1,0)[12], AIC=1377.52
(50) ARIMA(6,0,13)(0,1,0)[12], AIC=1377.52
(51) ARIMA(4,0,15)(0,1,0)[12], AIC=1377.61
(52) ARIMA(2,0,17)(0,1,0)[12], AIC=1377.67
(53) ARIMA(2,0,17)(0,1,0)[12], AIC=1377.67
(54) ARIMA(3,0,16)(0,1,0)[12], AIC=1377.71
(55) ARIMA(7,0,13)(0,1,0)[12], AIC=1378.15
(56) ARIMA(8,0,12)(0,1,0)[12], AIC=1378.47
(57) ARIMA(8,0,12)(0,1,0)[12], AIC=1378.47
(58) ARIMA(9,0,12)(0,1,0)[12], AIC=1378.79
(59) ARIMA(6,0,14)(0,1,0)[12], AIC=1378.97
(60) ARIMA(6,0,14)(0,1,0)[12], AIC=1378.97
(61) ARIMA(4,0,16)(0,1,0)[12], AIC=1379.46
(62) ARIMA(3,0,17)(0,1,0)[12], AIC=1379.51
(63) ARIMA(5,0,15)(0,1,0)[12], AIC=1379.68
(64) ARIMA(8,0,13)(0,1,0)[12], AIC=1380.15
(65) ARIMA(8,0,13)(0,1,0)[12], AIC=1380.15
(66) ARIMA(7,0,14)(0,1,0)[12], AIC=1380.53
(67) ARIMA(10,0,12)(0,1,0)[12], AIC=1380.62
(68) ARIMA(5,0,16)(0,1,0)[12], AIC=1380.90
(69) ARIMA(4,0,17)(0,1,0)[12], AIC=1381.04
(70) ARIMA(6,0,15)(0,1,0)[12], AIC=1381.49
(71) ARIMA(6,0,15)(0,1,0)[12], AIC=1381.49
(72) ARIMA(10,0,13)(0,1,0)[12], AIC=1381.66
(73) ARIMA(9,0,13)(0,1,0)[12], AIC=1381.75
(74) ARIMA(8,0,14)(0,1,0)[12], AIC=1382.31
(75) ARIMA(8,0,14)(0,1,0)[12], AIC=1382.31
(76) ARIMA(5,0,17)(0,1,0)[12], AIC=1382.35
(77) ARIMA(11,0,12)(0,1,0)[12], AIC=1382.48
(78) ARIMA(9,0,14)(0,1,0)[12], AIC=1382.58
(79) ARIMA(6,0,16)(0,1,0)[12], AIC=1382.64
(80) ARIMA(6,0,16)(0,1,0)[12], AIC=1382.64
(81) ARIMA(12,0,0)(0,1,1)[12], AIC=1382.66
(82) ARIMA(7,0,15)(0,1,0)[12], AIC=1382.73
(83) ARIMA(10,0,14)(0,1,0)[12], AIC=1383.53
(84) ARIMA(11,0,13)(0,1,0)[12], AIC=1383.82
(85) ARIMA(6,0,17)(0,1,0)[12], AIC=1384.31
(86) ARIMA(6,0,17)(0,1,0)[12], AIC=1384.31
(87) ARIMA(7,0,16)(0,1,0)[12], AIC=1384.33
(88) ARIMA(8,0,15)(0,1,0)[12], AIC=1384.52
(89) ARIMA(8,0,15)(0,1,0)[12], AIC=1384.52
(90) ARIMA(9,0,15)(0,1,0)[12], AIC=1384.53
(91) ARIMA(13,0,12)(0,1,0)[12], AIC=1384.56
(92) ARIMA(11,0,14)(0,1,0)[12], AIC=1384.73
(93) ARIMA(13,0,13)(0,1,0)[12], AIC=1385.12
(94) ARIMA(10,0,15)(0,1,0)[12], AIC=1385.45
(95) ARIMA(12,0,13)(0,1,0)[12], AIC=1385.52
(96) ARIMA(8,0,16)(0,1,0)[12], AIC=1385.88
(97) ARIMA(8,0,16)(0,1,0)[12], AIC=1385.88
(98) ARIMA(12,0,12)(0,1,0)[12], AIC=1386.15
(99) ARIMA(7,0,17)(0,1,0)[12], AIC=1386.37
(100) ARIMA(11,0,15)(0,1,0)[12], AIC=1386.42
(101) ARIMA(9,0,16)(0,1,0)[12], AIC=1386.43
(102) ARIMA(13,0,14)(0,1,0)[12], AIC=1386.89
(103) ARIMA(12,0,14)(0,1,0)[12], AIC=1386.90
(104) ARIMA(10,0,16)(0,1,0)[12], AIC=1387.11
(105) ARIMA(14,0,13)(0,1,0)[12], AIC=1387.36
(106) ARIMA(15,0,0)(0,1,1)[12], AIC=1387.68
(107) ARIMA(14,0,12)(0,1,0)[12], AIC=1387.85
(108) ARIMA(11,0,16)(0,1,0)[12], AIC=1387.97
(109) ARIMA(8,0,17)(0,1,0)[12], AIC=1388.02
(110) ARIMA(8,0,17)(0,1,0)[12], AIC=1388.02
(111) ARIMA(15,0,12)(0,1,0)[12], AIC=1388.42
(112) ARIMA(13,0,15)(0,1,0)[12], AIC=1388.76
(113) ARIMA(12,0,15)(0,1,0)[12], AIC=1388.84
(114) ARIMA(14,0,14)(0,1,0)[12], AIC=1388.89
(115) ARIMA(9,0,17)(0,1,0)[12], AIC=1388.95
(116) ARIMA(10,0,17)(0,1,0)[12], AIC=1389.04
(117) ARIMA(15,0,13)(0,1,0)[12], AIC=1389.30
(118) ARIMA(17,0,12)(0,1,0)[12], AIC=1389.57
(119) ARIMA(11,0,17)(0,1,0)[12], AIC=1389.72
(120) ARIMA(16,0,12)(0,1,0)[12], AIC=1389.87
(121) ARIMA(12,0,16)(0,1,0)[12], AIC=1390.33
(122) ARIMA(15,0,14)(0,1,0)[12], AIC=1390.73
(123) ARIMA(14,0,15)(0,1,0)[12], AIC=1390.78
(124) ARIMA(17,0,13)(0,1,0)[12], AIC=1390.81
(125) ARIMA(16,0,13)(0,1,0)[12], AIC=1390.83
(126) ARIMA(6,0,0)(0,1,1)[12], AIC=1391.61
(127) ARIMA(13,0,16)(0,1,0)[12], AIC=1391.95
(128) ARIMA(12,0,17)(0,1,0)[12], AIC=1392.26
(129) ARIMA(13,0,17)(0,1,0)[12], AIC=1392.28
(130) ARIMA(17,0,14)(0,1,0)[12], AIC=1392.49
(131) ARIMA(11,0,0)(0,1,4)[12], AIC=1392.63
(132) ARIMA(16,0,14)(0,1,0)[12], AIC=1392.69
(133) ARIMA(15,0,15)(0,1,0)[12], AIC=1392.85
(134) ARIMA(14,0,16)(0,1,0)[12], AIC=1392.97
(135) ARIMA(16,0,15)(0,1,0)[12], AIC=1394.72
(136) ARIMA(17,0,15)(0,1,0)[12], AIC=1394.84
(137) ARIMA(15,0,16)(0,1,0)[12], AIC=1395.10
(138) ARIMA(14,0,17)(0,1,0)[12], AIC=1395.38
(139) ARIMA(15,0,17)(0,1,0)[12], AIC=1396.22
(140) ARIMA(17,0,16)(0,1,0)[12], AIC=1396.75
(141) ARIMA(16,0,16)(0,1,0)[12], AIC=1397.15
(142) ARIMA(16,0,17)(0,1,0)[12], AIC=1398.28
(143) ARIMA(10,0,0)(0,1,4)[12], AIC=1398.90
(144) ARIMA(17,0,17)(0,1,0)[12], AIC=1399.30
(145) ARIMA(9,0,11)(0,1,0)[12], AIC=1405.22
(146) ARIMA(11,0,10)(0,1,0)[12], AIC=1405.63
(147) ARIMA(11,0,11)(0,1,0)[12], AIC=1406.36
(148) ARIMA(10,0,11)(0,1,0)[12], AIC=1411.83
(149) ARIMA(14,0,11)(0,1,0)[12], AIC=1413.61
(150) ARIMA(13,0,11)(0,1,0)[12], AIC=1415.15
(151) ARIMA(8,0,11)(0,1,0)[12], AIC=1415.83
(152) ARIMA(8,0,11)(0,1,0)[12], AIC=1415.83
(153) ARIMA(15,0,11)(0,1,0)[12], AIC=1418.52
(154) ARIMA(16,0,11)(0,1,0)[12], AIC=1418.96
(155) ARIMA(12,0,11)(0,1,0)[12], AIC=1419.07
(156) ARIMA(7,0,11)(0,1,0)[12], AIC=1420.16
(157) ARIMA(15,0,10)(0,1,0)[12], AIC=1427.03
(158) ARIMA(0,0,11)(3,1,0)[12], AIC=1427.35
(159) ARIMA(16,0,10)(0,1,0)[12], AIC=1428.10
(160) ARIMA(16,0,10)(0,1,0)[12], AIC=1428.10
(161) ARIMA(14,0,10)(0,1,0)[12], AIC=1429.08
(162) ARIMA(0,0,10)(3,1,0)[12], AIC=1430.51
(163) ARIMA(4,0,11)(0,1,0)[12], AIC=1430.74
(164) ARIMA(3,0,0)(0,1,1)[12], AIC=1431.92
(165) ARIMA(3,0,11)(0,1,0)[12], AIC=1435.07
(166) ARIMA(17,0,10)(0,1,0)[12], AIC=1436.99
(167) ARIMA(12,0,10)(0,1,0)[12], AIC=1437.18
(168) ARIMA(0,0,11)(2,1,0)[12], AIC=1443.12
(169) ARIMA(0,0,0)(0,1,3)[12], AIC=1443.50
(170) ARIMA(0,0,0)(0,1,1)[12], AIC=1443.59
(171) ARIMA(0,0,0)(0,1,1)[12], AIC=1443.59
(172) ARIMA(5,0,11)(0,1,0)[12], AIC=1443.88
(173) ARIMA(0,0,0)(2,1,1)[12], AIC=1444.11
(174) ARIMA(13,0,10)(0,1,0)[12], AIC=1444.13
(175) ARIMA(0,0,0)(1,1,2)[12], AIC=1445.15
(176) ARIMA(0,0,0)(1,1,3)[12], AIC=1445.28
(177) ARIMA(0,0,0)(0,1,4)[12], AIC=1445.49
(178) ARIMA(0,0,0)(0,1,2)[12], AIC=1445.59
(179) ARIMA(0,0,0)(1,1,1)[12], AIC=1445.59
(180) ARIMA(0,0,0)(4,1,3)[12], AIC=1445.70
(181) ARIMA(0,0,0)(4,1,4)[12], AIC=1445.74
(182) ARIMA(0,0,0)(2,1,2)[12], AIC=1445.92
(183) ARIMA(0,0,0)(4,1,2)[12], AIC=1445.95
(184) ARIMA(0,0,0)(3,1,1)[12], AIC=1446.07
(185) ARIMA(0,0,0)(2,1,3)[12], AIC=1446.38
(186) ARIMA(0,0,0)(1,1,4)[12], AIC=1446.57
(187) ARIMA(0,0,0)(4,1,1)[12], AIC=1446.70
(188) ARIMA(0,0,10)(2,1,0)[12], AIC=1446.97
(189) ARIMA(0,0,0)(3,1,2)[12], AIC=1446.99
(190) ARIMA(0,0,0)(3,1,3)[12], AIC=1447.78
(191) ARIMA(0,0,0)(2,1,4)[12], AIC=1447.80
(192) ARIMA(0,0,0)(4,1,0)[12], AIC=1448.24
(193) ARIMA(6,0,11)(0,1,0)[12], AIC=1449.08
(194) ARIMA(6,0,11)(0,1,0)[12], AIC=1449.08
(195) ARIMA(0,0,0)(3,1,4)[12], AIC=1449.21
(196) ARIMA(16,0,9)(0,1,0)[12], AIC=1450.39
(197) ARIMA(17,0,9)(0,1,0)[12], AIC=1450.46
(198) ARIMA(2,0,11)(0,1,0)[12], AIC=1451.97
(199) ARIMA(2,0,11)(0,1,0)[12], AIC=1451.97
(200) ARIMA(0,0,0)(3,1,0)[12], AIC=1454.13
(201) ARIMA(16,0,8)(0,1,0)[12], AIC=1460.41
(202) ARIMA(0,0,0)(2,1,0)[12], AIC=1464.14
(203) ARIMA(17,0,8)(0,1,0)[12], AIC=1464.61
(204) ARIMA(16,0,7)(0,1,0)[12], AIC=1469.14
(205) ARIMA(17,0,7)(0,1,0)[12], AIC=1469.22
(206) ARIMA(17,0,6)(0,1,0)[12], AIC=1471.02
(207) ARIMA(9,0,10)(0,1,0)[12], AIC=1471.14
(208) ARIMA(8,0,10)(0,1,0)[12], AIC=1478.49
(209) ARIMA(8,0,10)(0,1,0)[12], AIC=1478.49
(210) ARIMA(16,0,5)(0,1,0)[12], AIC=1479.55
(211) ARIMA(16,0,6)(0,1,0)[12], AIC=1480.98
(212) ARIMA(10,0,10)(0,1,0)[12], AIC=1495.41
(213) ARIMA(17,0,5)(0,1,0)[12], AIC=1496.61
(214) ARIMA(6,0,10)(0,1,0)[12], AIC=1497.53
(215) ARIMA(6,0,10)(0,1,0)[12], AIC=1497.53
(216) ARIMA(1,0,11)(0,1,0)[12], AIC=1503.08
(217) ARIMA(0,0,11)(1,1,0)[12], AIC=1511.21
(218) ARIMA(0,0,6)(1,1,0)[12], AIC=1514.47
(219) ARIMA(0,0,8)(1,1,0)[12], AIC=1516.15
(220) ARIMA(0,0,10)(1,1,0)[12], AIC=1516.98
(221) ARIMA(12,0,0)(0,1,0)[12], AIC=1522.91
(222) ARIMA(0,0,0)(1,1,0)[12], AIC=1526.64
(223) ARIMA(5,0,10)(0,1,0)[12], AIC=1526.92
(224) ARIMA(10,0,6)(0,1,0)[12], AIC=1541.36
(225) ARIMA(4,0,10)(0,1,0)[12], AIC=1542.28
(226) ARIMA(7,0,10)(0,1,0)[12], AIC=1546.44
(227) ARIMA(3,0,10)(0,1,0)[12], AIC=1552.06
(228) ARIMA(2,0,10)(0,1,0)[12], AIC=1552.72
(229) ARIMA(2,0,10)(0,1,0)[12], AIC=1552.72
(230) ARIMA(2,0,10)(0,1,0)[12], AIC=1552.72
(231) ARIMA(2,0,9)(0,1,0)[12], AIC=1555.97
(232) ARIMA(2,0,8)(0,1,0)[12], AIC=1559.09
(233) ARIMA(1,0,8)(0,1,0)[12], AIC=1593.01
(234) ARIMA(0,0,11)(0,1,0)[12], AIC=1599.58
(235) ARIMA(2,0,5)(0,1,0)[12], AIC=1601.11
(236) ARIMA(2,0,6)(0,1,0)[12], AIC=1621.74
(237) ARIMA(0,0,9)(0,1,0)[12], AIC=1621.79
(238) ARIMA(0,0,10)(0,1,0)[12], AIC=1622.73
(239) ARIMA(0,0,10)(0,1,0)[12], AIC=1622.73
(240) ARIMA(1,0,9)(0,1,0)[12], AIC=1623.45
(241) ARIMA(1,0,10)(0,1,0)[12], AIC=1624.20
(242) ARIMA(1,0,10)(0,1,0)[12], AIC=1624.20
(243) ARIMA(1,0,5)(0,1,0)[12], AIC=1627.56
(244) ARIMA(1,0,6)(0,1,0)[12], AIC=1630.49
(245) ARIMA(0,0,7)(0,1,0)[12], AIC=1630.59
(246) ARIMA(0,0,8)(0,1,0)[12], AIC=1630.74
(247) ARIMA(1,0,7)(0,1,0)[12], AIC=1631.92
(248) ARIMA(2,0,7)(0,1,0)[12], AIC=1633.08
(249) ARIMA(0,0,6)(0,1,0)[12], AIC=1635.97
(250) ARIMA(6,0,0)(0,1,0)[12], AIC=1650.81
(251) ARIMA(0,0,5)(0,1,0)[12], AIC=1651.91
(252) ARIMA(0,0,0)(0,1,0)[12], AIC=1653.67

MODELOS PROPUESTOS¶

In [89]:
from statsmodels.tsa.statespace.sarimax import SARIMAX
from verificacion_de_supuestos import * # Una funcion que hice para verificar si los parametros son estacionarios/invertibles
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import ttest_1samp
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.tools.tools import add_constant
from scipy.stats import jarque_bera
from statsmodels.stats.diagnostic import lilliefors

MODELOS PARA LA SERIE SIN TRANSFORMAR¶

1 MODELO¶

In [36]:
modelo=SARIMAX(zpre,
               order=(12,0,0),
               seasonal_order=(0,1,4,12)).fit()

modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
Out[36]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(12, 0, 0)x(0, 1, [1, 2, 3, 4], 12) Log Likelihood -666.522
Date: Sun, 27 Apr 2025 AIC 1367.044
Time: 04:48:12 BIC 1443.100
Sample: 0 HQIC 1396.549
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.1079 0.032 3.338 0.001 0.045 0.171
ar.L2 -0.0271 0.037 -0.738 0.460 -0.099 0.045
ar.L3 0.0504 0.037 1.367 0.172 -0.022 0.123
ar.L4 -0.0689 0.042 -1.634 0.102 -0.152 0.014
ar.L5 -0.1086 0.050 -2.193 0.028 -0.206 -0.012
ar.L6 -0.1103 0.059 -1.872 0.061 -0.226 0.005
ar.L7 -0.0049 0.050 -0.097 0.923 -0.104 0.094
ar.L8 -0.0450 0.043 -1.039 0.299 -0.130 0.040
ar.L9 0.0411 0.035 1.183 0.237 -0.027 0.109
ar.L10 -0.0133 0.033 -0.406 0.685 -0.078 0.051
ar.L11 0.0968 0.036 2.662 0.008 0.026 0.168
ar.L12 0.5442 0.074 7.328 0.000 0.399 0.690
ma.S.L12 -1.3768 0.084 -16.419 0.000 -1.541 -1.212
ma.S.L24 0.3656 0.079 4.639 0.000 0.211 0.520
ma.S.L36 0.1054 0.060 1.756 0.079 -0.012 0.223
ma.S.L48 -0.0790 0.039 -2.042 0.041 -0.155 -0.003
sigma2 0.4378 0.024 18.022 0.000 0.390 0.485
Ljung-Box (L1) (Q): 0.34 Jarque-Bera (JB): 376.73
Prob(Q): 0.56 Prob(JB): 0.00
Heteroskedasticity (H): 0.82 Skew: 0.89
Prob(H) (two-sided): 0.14 Kurtosis: 6.28


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

PRINCIPIO DE PARSIMONIA¶

In [37]:
parsimonia(modelo) 
Hay coeficientes no significativos, no se cumple el principio de parsimonia

MODELO ADMISIBLE¶

In [38]:
estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.10787169558218329, 'p2': 0.027085643770389415, 'p3': -0.050379964730704004, 'p4': 0.06894551834156623, 'p5': 0.10861700542171245, 'p6': 0.11029529000485752, 'p7': 0.00488096912398539, 'p8': 0.04497037748451452, 'p9': -0.04112819961848579, 'p10': 0.013332216884030874, 'p11': -0.09678291912170793, 'p12': -0.5442160802231802}

El grado del polinomio es: 12

Raíces del polinomio característico: [-1.10915274+0.j         -0.92789955+0.52544644j -0.92789955-0.52544644j
 -0.52132244+0.91355779j -0.52132244-0.91355779j -0.01548825+1.04638174j
 -0.01548825-1.04638174j  0.53546751+0.91668064j  0.53546751-0.91668064j
  1.05996383+0.j          0.86491762+0.50841197j  0.86491762-0.50841197j]

Módulo de las raíces: [1.10915274 1.06634494 1.06634494 1.05183883 1.05183883 1.04649636
 1.04649636 1.06161615 1.06161615 1.05996383 1.00327724 1.00327724]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)
In [39]:
invertible(modelo)
El polinomio de la parte media móvil es: {'p12': -1.376803188196061, 'p24': 0.36563094250781014, 'p36': 0.10536384747112354, 'p48': -0.07904343456059729}

El grado del polinomio es: 48

Raíces del polinomio característico: [-1.05513205e+00+0.28272178j -1.05513205e+00-0.28272178j
 -1.06072484e+00+0.06092667j -1.06072484e+00-0.06092667j
 -1.00194881e+00+0.j         -7.72410271e-01+0.77241027j
 -7.72410271e-01-0.77241027j -9.49077989e-01+0.47759838j
 -9.49077989e-01-0.47759838j -8.88151320e-01+0.58312646j
 -8.88151320e-01-0.58312646j -8.67713125e-01+0.50097441j
 -8.67713125e-01-0.50097441j -5.83126461e-01+0.88815132j
 -5.83126461e-01-0.88815132j -4.77598375e-01+0.94907799j
 -4.77598375e-01-0.94907799j -5.00974406e-01+0.86771313j
 -5.00974406e-01-0.86771313j -2.82721781e-01+1.05513205j
 -2.82721781e-01-1.05513205j -6.09266684e-02+1.06072484j
 -6.09266684e-02-1.06072484j  6.09266684e-02+1.06072484j
  6.09266684e-02-1.06072484j  4.53730117e-17+1.00194881j
  4.53730117e-17-1.00194881j  2.82721781e-01+1.05513205j
  2.82721781e-01-1.05513205j  4.77598375e-01+0.94907799j
  4.77598375e-01-0.94907799j  5.83126461e-01+0.88815132j
  5.83126461e-01-0.88815132j  5.00974406e-01+0.86771313j
  5.00974406e-01-0.86771313j  7.72410271e-01+0.77241027j
  7.72410271e-01-0.77241027j  1.05513205e+00+0.28272178j
  1.05513205e+00-0.28272178j  1.06072484e+00+0.06092667j
  1.06072484e+00-0.06092667j  1.00194881e+00+0.j
  8.88151320e-01+0.58312646j  8.88151320e-01-0.58312646j
  9.49077989e-01+0.47759838j  9.49077989e-01-0.47759838j
  8.67713125e-01+0.50097441j  8.67713125e-01-0.50097441j]

Módulo de las raíces: [1.09235308 1.09235308 1.06247317 1.06247317 1.00194881 1.09235308
 1.09235308 1.06247317 1.06247317 1.06247317 1.06247317 1.00194881
 1.00194881 1.06247317 1.06247317 1.06247317 1.06247317 1.00194881
 1.00194881 1.09235308 1.09235308 1.06247317 1.06247317 1.06247317
 1.06247317 1.00194881 1.00194881 1.09235308 1.09235308 1.06247317
 1.06247317 1.06247317 1.06247317 1.00194881 1.00194881 1.09235308
 1.09235308 1.09235308 1.09235308 1.06247317 1.06247317 1.00194881
 1.06247317 1.06247317 1.06247317 1.06247317 1.00194881 1.00194881]

¿Las raíces están fuera del círculo unitario?  True

El modelo es invertible

:)

RESIDUOS INDEPENDIENTES¶

La prueba de Ljung-Box evalúa si los residuos de un modelo están autocorrelacionados, es decir, si existe dependencia temporal que el modelo no ha capturado. El resumen muestra lo siguiente:

Estadístico Q de Ljung-Box (L=1):  0.00  
Valor-p asociado (Prob(Q)):       0.95
  • El estadístico Q = 0.00 indica que no se detecta autocorrelación en los residuos al primer rezago.
  • El valor-p = 0.95 es mucho mayor que el umbral común de significancia (α = 0.05).

Por lo tanto, no se rechaza la hipótesis nula de independencia, los residuos se comportan de manera independiente.

La prueba Ljung-Box evalúa la autocorrelación de los residuos, cuando estos residuos son residuales de un modelo ajustado, el estadístico debe corregirse para reflejar la pérdida de grados de libertad debido a la estimación de parámetros.

La prueba original de Ljung-Box asume que los datos son observaciones independientes sin parámetros estimados Cuando se usan residuos ajustados por un modelo con $k$ parámetros estimados, los grados de libertad del estadístico Q se reducen en función de estos parámetros. El argumento model_df en acorr_ljungbox() resta esos grados de libertad automáticamente.

In [40]:
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
Out[40]:
lb_stat lb_pvalue
1 1.410877 NaN
2 1.666736 NaN
3 1.668529 NaN
4 2.153928 NaN
5 3.756515 NaN
6 4.271537 NaN
7 4.271537 NaN
8 4.398558 NaN
9 4.673186 NaN
10 5.854609 NaN
11 6.235510 NaN
12 6.364705 NaN
13 8.908873 NaN
14 9.013624 NaN
15 9.115022 NaN
16 9.301973 NaN
17 9.622827 0.001922
18 9.634894 0.008087
19 9.648719 0.021801
20 9.650120 0.046752
21 10.177061 0.070371
22 10.403762 0.108646
23 10.403778 0.166823
24 10.406790 0.237627
25 10.406958 0.318555
26 10.956519 0.360916
27 11.577235 0.396244
28 11.768293 0.464463
29 12.219944 0.509696
30 12.223270 0.588378
31 12.229827 0.661556
32 14.130904 0.588962
33 14.189628 0.653636
34 16.160520 0.581346
35 16.186598 0.644794
36 16.663770 0.674686
37 16.699639 0.729145
38 17.046283 0.760802
39 17.392114 0.789577
40 18.027704 0.801662
41 19.922150 0.750897
42 20.004090 0.791363
43 20.038756 0.829126
44 20.096537 0.860920
45 20.910376 0.862394
46 21.002354 0.887814
47 21.052280 0.910525
48 21.957194 0.908534

Hay dependencia en los primeros dos residuos.

RESIDUOS CON MEDIA CERO¶

In [41]:
modelo.resid.mean()
Out[41]:
-0.04082317526120878
In [42]:
ttest_1samp(modelo.resid, 0)
Out[42]:
TtestResult(statistic=-1.5347957438291504, pvalue=0.12531394692788528, df=659)

La media de los residuos puede ser 0.

RESIDUOS CON VARIANZA CONSTANTE¶

In [43]:
het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))
Out[43]:
(2.3754715860699527,
 0.12325437239232831,
 2.376827864684162,
 0.12362813746633024)

Los residuos tienen varianza constante.

RESIDUOS CON DISTRIBUCIÓN NORMAL¶

In [44]:
print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))
SignificanceResult(statistic=417.01556109393744, pvalue=2.793968301284977e-91)
(0.11546272377541411, 0.0009999999999998899)

Los residuos no siguen una distribución normal.

In [45]:
normal_relajacion(modelo.resid)
75.30% de los residuos están dentro de ±1σ (esperado ≈ 68%)
93.94% de los residuos están dentro de ±2σ (esperado ≈ 95%)
98.79% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)

Pero son bastante parecidos.

GRÁFICO DE RESIDUOS¶

In [46]:
from matplotlib.collections import PathCollection
In [47]:
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 4, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-01-analisis-de-residuos-m1.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

2 MODELO¶

In [48]:
modelo=SARIMAX(zpre,
               order=(13,0,0),
               seasonal_order=(0,1,4,12)).fit()

modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
Out[48]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(13, 0, 0)x(0, 1, [1, 2, 3, 4], 12) Log Likelihood -665.685
Date: Sun, 27 Apr 2025 AIC 1367.369
Time: 04:49:22 BIC 1447.899
Sample: 0 HQIC 1398.609
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.1369 0.035 3.963 0.000 0.069 0.205
ar.L2 -0.0190 0.036 -0.522 0.602 -0.090 0.052
ar.L3 0.0501 0.037 1.363 0.173 -0.022 0.122
ar.L4 -0.0697 0.042 -1.662 0.097 -0.152 0.013
ar.L5 -0.1117 0.050 -2.254 0.024 -0.209 -0.015
ar.L6 -0.1101 0.060 -1.839 0.066 -0.227 0.007
ar.L7 -0.0104 0.051 -0.205 0.838 -0.110 0.089
ar.L8 -0.0511 0.044 -1.170 0.242 -0.137 0.034
ar.L9 0.0331 0.035 0.944 0.345 -0.036 0.102
ar.L10 -0.0093 0.033 -0.281 0.778 -0.074 0.055
ar.L11 0.1034 0.037 2.802 0.005 0.031 0.176
ar.L12 0.5385 0.080 6.699 0.000 0.381 0.696
ar.L13 -0.0460 0.041 -1.115 0.265 -0.127 0.035
ma.S.L12 -1.3756 0.100 -13.725 0.000 -1.572 -1.179
ma.S.L24 0.3674 0.086 4.275 0.000 0.199 0.536
ma.S.L36 0.0973 0.060 1.635 0.102 -0.019 0.214
ma.S.L48 -0.0784 0.039 -1.994 0.046 -0.155 -0.001
sigma2 0.4340 0.030 14.503 0.000 0.375 0.493
Ljung-Box (L1) (Q): 0.04 Jarque-Bera (JB): 347.08
Prob(Q): 0.85 Prob(JB): 0.00
Heteroskedasticity (H): 0.83 Skew: 0.87
Prob(H) (two-sided): 0.17 Kurtosis: 6.14


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

PRINCIPIO DE PARSIMONIA¶

In [49]:
parsimonia(modelo) 
Hay coeficientes no significativos, no se cumple el principio de parsimonia

MODELO ADMISIBLE¶

In [50]:
estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.13691650540076578, 'p2': 0.01901377784903687, 'p3': -0.0500667165465603, 'p4': 0.069698229094628, 'p5': 0.11171395157927178, 'p6': 0.11009957522513744, 'p7': 0.010425775539985482, 'p8': 0.051137304909435626, 'p9': -0.03308758191153757, 'p10': 0.009267663653695928, 'p11': -0.1033592404899208, 'p12': -0.5385061478654415, 'p13': 0.046023546700947306}

El grado del polinomio es: 13

Raíces del polinomio característico: [11.88851581+0.j         -1.09964507+0.j         -0.92538665+0.52179453j
 -0.92538665-0.52179453j -0.52533631+0.91031656j -0.52533631-0.91031656j
 -0.02218345+1.04735764j -0.02218345-1.04735764j  0.5303036 +0.92193348j
  0.5303036 -0.92193348j  1.06751362+0.j          0.86474363+0.50810694j
  0.86474363-0.50810694j]

Módulo de las raíces: [11.88851581  1.09964507  1.06236057  1.06236057  1.05102544  1.05102544
  1.04759254  1.04759254  1.06357099  1.06357099  1.06751362  1.00297268
  1.00297268]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)
In [51]:
invertible(modelo)
El polinomio de la parte media móvil es: {'p12': -1.3756005193776482, 'p24': 0.3674412910433979, 'p36': 0.09734997881875754, 'p48': -0.07835471066657075}

El grado del polinomio es: 48

Raíces del polinomio característico: [-1.05628158e+00+0.2830298j  -1.05628158e+00-0.2830298j
 -1.06074970e+00+0.06258621j -1.06074970e+00-0.06258621j
 -1.00135998e+00+0.j         -7.73251787e-01+0.77325179j
 -7.73251787e-01-0.77325179j -9.49929287e-01+0.4761736j
 -9.49929287e-01-0.4761736j  -8.87343081e-01+0.58457609j
 -8.87343081e-01-0.58457609j -8.67203179e-01+0.50067999j
 -8.67203179e-01-0.50067999j -5.84576092e-01+0.88734308j
 -5.84576092e-01-0.88734308j -4.76173604e-01+0.94992929j
 -4.76173604e-01-0.94992929j -5.00679989e-01+0.86720318j
 -5.00679989e-01-0.86720318j -2.83029797e-01+1.05628158j
 -2.83029797e-01-1.05628158j -6.25862060e-02+1.0607497j
 -6.25862060e-02-1.0607497j   6.25862060e-02+1.0607497j
  6.25862060e-02-1.0607497j   2.16421979e-16+1.00135998j
  2.16421979e-16-1.00135998j  2.83029797e-01+1.05628158j
  2.83029797e-01-1.05628158j  4.76173604e-01+0.94992929j
  4.76173604e-01-0.94992929j  5.84576092e-01+0.88734308j
  5.84576092e-01-0.88734308j  5.00679989e-01+0.86720318j
  5.00679989e-01-0.86720318j  7.73251787e-01+0.77325179j
  7.73251787e-01-0.77325179j  1.05628158e+00+0.2830298j
  1.05628158e+00-0.2830298j   1.06074970e+00+0.06258621j
  1.06074970e+00-0.06258621j  1.00135998e+00+0.j
  8.87343081e-01+0.58457609j  8.87343081e-01-0.58457609j
  9.49929287e-01+0.4761736j   9.49929287e-01-0.4761736j
  8.67203179e-01+0.50067999j  8.67203179e-01-0.50067999j]

Módulo de las raíces: [1.09354316 1.09354316 1.06259444 1.06259444 1.00135998 1.09354316
 1.09354316 1.06259444 1.06259444 1.06259444 1.06259444 1.00135998
 1.00135998 1.06259444 1.06259444 1.06259444 1.06259444 1.00135998
 1.00135998 1.09354316 1.09354316 1.06259444 1.06259444 1.06259444
 1.06259444 1.00135998 1.00135998 1.09354316 1.09354316 1.06259444
 1.06259444 1.06259444 1.06259444 1.00135998 1.00135998 1.09354316
 1.09354316 1.09354316 1.09354316 1.06259444 1.06259444 1.00135998
 1.06259444 1.06259444 1.06259444 1.06259444 1.00135998 1.00135998]

¿Las raíces están fuera del círculo unitario?  True

El modelo es invertible

:)

RESIDUOS INDEPENDIENTES¶

In [52]:
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
Out[52]:
lb_stat lb_pvalue
1 0.169555 NaN
2 0.204491 NaN
3 0.207502 NaN
4 0.727124 NaN
5 2.068440 NaN
6 2.452464 NaN
7 2.514097 NaN
8 2.809591 NaN
9 3.353825 NaN
10 4.344480 NaN
11 4.491955 NaN
12 4.756059 NaN
13 5.816217 NaN
14 6.013503 NaN
15 6.142261 NaN
16 6.454918 NaN
17 6.774386 NaN
18 6.783618 0.009200
19 6.786944 0.033592
20 6.809256 0.078232
21 7.196463 0.125863
22 7.428032 0.190703
23 7.441854 0.281909
24 7.447962 0.383773
25 7.506681 0.483077
26 8.263074 0.507868
27 8.898238 0.541788
28 9.010612 0.620913
29 9.533962 0.656772
30 9.544856 0.730662
31 9.545139 0.794626
32 11.177851 0.739888
33 11.194537 0.797321
34 13.243411 0.719745
35 13.251323 0.776431
36 13.698894 0.800953
37 13.804935 0.840247
38 14.263347 0.858024
39 14.567144 0.880135
40 15.189601 0.887710
41 16.919636 0.852069
42 16.949758 0.883608
43 16.977019 0.909776
44 17.083844 0.929034
45 17.970345 0.926894
46 18.039168 0.943470
47 18.078955 0.957241
48 18.964369 0.955485

Hay dependencia en los primeros dos residuos.

RESIDUOS CON MEDIA CERO¶

In [53]:
modelo.resid.mean()
Out[53]:
-0.04384260119451572
In [54]:
ttest_1samp(modelo.resid, 0)
Out[54]:
TtestResult(statistic=-1.6518512196908124, pvalue=0.09904125403907901, df=659)

La media de los residuos puede ser 0.

RESIDUOS CON VARIANZA CONSTANTE¶

In [55]:
het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))
Out[55]:
(2.2724202057175913,
 0.1316941016400431,
 2.273361405689974,
 0.13209411797812018)

Los residuos tienen varianza constante.

RESIDUOS CON DISTRIBUCIÓN NORMAL¶

In [56]:
print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))
SignificanceResult(statistic=380.3670597932282, pvalue=2.537128226453911e-83)
(0.11315000343530657, 0.0009999999999998899)

Los residuos no siguen una distribución normal.

In [57]:
normal_relajacion(modelo.resid)
75.15% de los residuos están dentro de ±1σ (esperado ≈ 68%)
93.79% de los residuos están dentro de ±2σ (esperado ≈ 95%)
98.94% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)

Pero son bastante parecidos.

GRÁFICO DE RESIDUOS¶

In [58]:
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    # line.set_markerfacecolor('green')
    # line.set_markeredgecolor('white')
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(13, 0, 0)(0, 1, 4, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-02-analisis-de-residuos-m2.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

3 MODELO¶

In [59]:
modelo=SARIMAX(zpre,
               order=(12,0,0),
               seasonal_order=(0,1,2,12)).fit()

modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
Out[59]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(12, 0, 0)x(0, 1, [1, 2], 12) Log Likelihood -668.499
Date: Sun, 27 Apr 2025 AIC 1366.998
Time: 04:49:32 BIC 1434.107
Sample: 0 HQIC 1393.032
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.1241 0.033 3.797 0.000 0.060 0.188
ar.L2 -0.0278 0.037 -0.746 0.455 -0.101 0.045
ar.L3 0.0593 0.036 1.629 0.103 -0.012 0.131
ar.L4 -0.0777 0.043 -1.819 0.069 -0.161 0.006
ar.L5 -0.1077 0.050 -2.137 0.033 -0.206 -0.009
ar.L6 -0.1154 0.060 -1.914 0.056 -0.234 0.003
ar.L7 -0.0107 0.050 -0.214 0.831 -0.108 0.087
ar.L8 -0.0498 0.044 -1.129 0.259 -0.136 0.037
ar.L9 0.0450 0.036 1.263 0.207 -0.025 0.115
ar.L10 -0.0087 0.033 -0.263 0.793 -0.074 0.056
ar.L11 0.1075 0.035 3.042 0.002 0.038 0.177
ar.L12 0.5090 0.073 6.931 0.000 0.365 0.653
ma.S.L12 -1.3374 0.082 -16.355 0.000 -1.498 -1.177
ma.S.L24 0.3550 0.075 4.760 0.000 0.209 0.501
sigma2 0.4394 0.022 20.157 0.000 0.397 0.482
Ljung-Box (L1) (Q): 0.11 Jarque-Bera (JB): 385.54
Prob(Q): 0.74 Prob(JB): 0.00
Heteroskedasticity (H): 0.82 Skew: 0.91
Prob(H) (two-sided): 0.15 Kurtosis: 6.31


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

PRINCIPIO DE PARSIMONIA¶

In [60]:
parsimonia(modelo) 
Hay coeficientes no significativos, no se cumple el principio de parsimonia

MODELO ADMISIBLE¶

In [61]:
estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.12409609720857327, 'p2': 0.02784753160299208, 'p3': -0.05931370060464292, 'p4': 0.07773595850485018, 'p5': 0.10769845417883445, 'p6': 0.1153846887487895, 'p7': 0.010651434003279378, 'p8': 0.04983749921571501, 'p9': -0.04496683215549213, 'p10': 0.008718941011892536, 'p11': -0.10747007955137028, 'p12': -0.5090305660410539}

El grado del polinomio es: 12

Raíces del polinomio característico: [-1.12406749+0.j         -0.93496549+0.52952748j -0.93496549-0.52952748j
 -0.52508333+0.91815394j -0.52508333-0.91815394j -0.01640444+1.05341344j
 -0.01640444-1.05341344j  0.53677332+0.92317355j  0.53677332-0.92317355j
  1.06348079+0.j          0.8644098 +0.50791775j  0.8644098 -0.50791775j]

Módulo de las raíces: [1.12406749 1.07450445 1.07450445 1.05769521 1.05769521 1.05354116
 1.05354116 1.06788342 1.06788342 1.06348079 1.00258902 1.00258902]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)
In [62]:
invertible(modelo)
El polinomio de la parte media móvil es: {'p12': -1.3373896067612656, 'p24': 0.3549819568174022}

El grado del polinomio es: 24

Raíces del polinomio característico: [-1.08759159e+00+0.j         -1.00234440e+00+0.j
 -9.41881945e-01+0.54379579j -9.41881945e-01-0.54379579j
 -8.68055716e-01+0.5011722j  -8.68055716e-01-0.5011722j
 -5.43795795e-01+0.94188195j -5.43795795e-01-0.94188195j
 -5.01172202e-01+0.86805572j -5.01172202e-01-0.86805572j
 -7.34340403e-16+1.08759159j -7.34340403e-16-1.08759159j
  9.15927823e-17+1.0023444j   9.15927823e-17-1.0023444j
  5.43795795e-01+0.94188195j  5.43795795e-01-0.94188195j
  5.01172202e-01+0.86805572j  5.01172202e-01-0.86805572j
  9.41881945e-01+0.54379579j  9.41881945e-01-0.54379579j
  1.08759159e+00+0.j          1.00234440e+00+0.j
  8.68055716e-01+0.5011722j   8.68055716e-01-0.5011722j ]

Módulo de las raíces: [1.08759159 1.0023444  1.08759159 1.08759159 1.0023444  1.0023444
 1.08759159 1.08759159 1.0023444  1.0023444  1.08759159 1.08759159
 1.0023444  1.0023444  1.08759159 1.08759159 1.0023444  1.0023444
 1.08759159 1.08759159 1.08759159 1.0023444  1.0023444  1.0023444 ]

¿Las raíces están fuera del círculo unitario?  True

El modelo es invertible

:)

RESIDUOS INDEPENDIENTES¶

In [63]:
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
Out[63]:
lb_stat lb_pvalue
1 0.852590 NaN
2 1.026959 NaN
3 1.042580 NaN
4 1.330274 NaN
5 2.702490 NaN
6 2.993605 NaN
7 3.102314 NaN
8 3.441180 NaN
9 3.616131 NaN
10 4.475840 NaN
11 4.662731 NaN
12 4.990946 NaN
13 7.916099 NaN
14 8.223545 NaN
15 8.469831 0.003611
16 8.751457 0.012579
17 9.104454 0.027934
18 9.135493 0.057801
19 9.135514 0.103780
20 9.139624 0.165877
21 9.745116 0.203483
22 10.168911 0.253370
23 10.181542 0.335989
24 10.580242 0.391141
25 10.587480 0.478441
26 11.470408 0.489092
27 12.051651 0.523413
28 12.201265 0.590143
29 12.738996 0.622450
30 12.744188 0.691362
31 12.749628 0.752779
32 14.953637 0.665148
33 14.972683 0.724332
34 17.346330 0.630384
35 17.349312 0.689726
36 18.361216 0.684370
37 18.361346 0.737605
38 18.844803 0.760215
39 19.320376 0.781506
40 19.893036 0.796598
41 21.815050 0.746715
42 21.912901 0.785308
43 21.985648 0.820783
44 22.126466 0.849404
45 23.141840 0.843941
46 23.250141 0.870261
47 23.337109 0.893652
48 23.478694 0.912247

Hay dependencia en los primeros dos residuos.

RESIDUOS CON MEDIA CERO¶

In [64]:
modelo.resid.mean()
Out[64]:
-0.03975299918077584
In [65]:
ttest_1samp(modelo.resid, 0)
Out[65]:
TtestResult(statistic=-1.4907248613412258, pvalue=0.1365122225435534, df=659)

La media de los residuos puede ser 0.

RESIDUOS CON VARIANZA CONSTANTE¶

In [66]:
het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))
Out[66]:
(2.3477675031968714,
 0.12546257027135715,
 2.3490090062319307,
 0.1258433113034253)

Los residuos tienen varianza constante.

RESIDUOS CON DISTRIBUCIÓN NORMAL¶

In [67]:
print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))
SignificanceResult(statistic=421.16509310160836, pvalue=3.508827740565011e-92)
(0.10731720679113188, 0.0009999999999998899)

Los residuos no siguen una distribución normal.

In [68]:
normal_relajacion(modelo.resid)
75.15% de los residuos están dentro de ±1σ (esperado ≈ 68%)
94.09% de los residuos están dentro de ±2σ (esperado ≈ 95%)
98.79% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)

Pero son bastante parecidos.

GRÁFICO DE RESIDUOS¶

In [69]:
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    # line.set_markerfacecolor('green')
    # line.set_markeredgecolor('white')
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

OTROS MODELOS QUE SE CREYÓ QUE PODÍAN SER¶

Luego de ver la FAC y FACP:

In [70]:
fac = FAC(len(dzpre), acf(dzpre, nlags=12*12)[1:] )
Valores de autocorrelacion significativos:
r1: 0.08495427799538309
r5: -0.09455404757011152
r12: -0.41418332950296094
r69: 0.13811523521698316
r73: -0.11256100076387932
r82: -0.09936326169869977
r139: 0.10081286249593209
r144: -0.103722152043645
In [71]:
facp = FACP(len(dzpre), pacf(dzpre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos:
rho 1: 0.08508558290727704
rho 5: -0.08389897753926691
rho 12: -0.43583220424270375
rho 17: -0.12023508334631484
rho 24: -0.30097443805877616
rho 29: -0.10198920198664936
rho 36: -0.15950168604080825
rho 48: -0.14036369491687153
rho 57: -0.09168717760514207
rho 60: -0.1239606605892366
rho 61: 0.13800141862943802
rho 69: 0.10269095820413282
rho 72: -0.09462958643919699
rho 82: -0.08002594037040578
rho 85: -0.08198325512480845
rho 93: 0.092203788733838
rho 109: 0.09582379018763891
rho 113: -0.09105599732819755
rho 119: 0.09171771966671013
rho 120: -0.13558724808031059

Se pensó que podía ser este modelo, pero es muy malo.

In [72]:
modelo=SARIMAX(zpre,
               order=([1,5],0,[1,5]),
               seasonal_order=([1],1,[1,2],12)).fit()

modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
  warn('Non-stationary starting autoregressive parameters'
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
  warn('Non-invertible starting MA parameters found.'
Out[72]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX([1, 5], 0, [1, 5])x(1, 1, 2, 12) Log Likelihood -5.099
Date: Sun, 27 Apr 2025 AIC 26.198
Time: 04:51:06 BIC 61.990
Sample: 0 HQIC 40.083
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 1.1476 32.468 0.035 0.972 -62.489 64.784
ar.L5 -0.8775 38.486 -0.023 0.982 -76.308 74.553
ma.L1 -0.9317 23.880 -0.039 0.969 -47.736 45.872
ma.L5 0.3314 36.742 0.009 0.993 -71.681 72.344
ar.S.L12 -0.3697 133.649 -0.003 0.998 -262.317 261.577
ma.S.L12 -0.1630 41.007 -0.004 0.997 -80.536 80.210
ma.S.L24 -0.7558 42.337 -0.018 0.986 -83.735 82.223
sigma2 0.4434 8.359 0.053 0.958 -15.941 16.827
Ljung-Box (L1) (Q): 294.50 Jarque-Bera (JB): 2134861.83
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.00 Skew: 15.72
Prob(H) (two-sided): 0.00 Kurtosis: 282.43


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 6.39e+17. Standard errors may be unstable.
In [73]:
modelo.resid.mean()
Out[73]:
-7.786026841214317e+48

También este, pero no es invertible

In [74]:
modelo=SARIMAX(zpre,
               order=(0,0, [5]),
               seasonal_order=(1,1,[1,2],12)).fit()

modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
Out[74]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(0, 0, [5])x(1, 1, [1, 2], 12) Log Likelihood -708.534
Date: Sun, 27 Apr 2025 AIC 1427.068
Time: 04:51:13 BIC 1449.437
Sample: 0 HQIC 1435.745
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ma.L5 -0.1707 0.044 -3.861 0.000 -0.257 -0.084
ar.S.L12 -0.8858 0.147 -6.011 0.000 -1.175 -0.597
ma.S.L12 0.2236 0.140 1.599 0.110 -0.050 0.498
ma.S.L24 -0.6343 0.088 -7.209 0.000 -0.807 -0.462
sigma2 0.5151 0.019 27.358 0.000 0.478 0.552
Ljung-Box (L1) (Q): 10.39 Jarque-Bera (JB): 283.30
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.83 Skew: 0.66
Prob(H) (two-sided): 0.18 Kurtosis: 5.96


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
In [75]:
invertible(modelo)
El polinomio de la parte media móvil es: {'p5': -0.17066687867489047, 'p12': 0.2235561375559531, 'p24': -0.6342522946355272}

El grado del polinomio es: 24

Raíces del polinomio característico: [-1.03809773+0.j         -0.97339645+0.26906311j -0.97339645-0.26906311j
 -0.88525875+0.51576505j -0.88525875-0.51576505j -0.71259807+0.7035491j
 -0.71259807-0.7035491j  -0.52270681+0.89327539j -0.52270681-0.89327539j
 -0.26081491+0.98110341j -0.26081491-0.98110341j  0.00733524+1.03148712j
  0.00733524-1.03148712j  1.02333621+0.j          0.97341295+0.25216589j
  0.97341295-0.25216589j  0.89999742+0.51572558j  0.89999742-0.51572558j
  0.71256537+0.7204749j   0.71256537-0.7204749j   0.50801366+0.89331475j
  0.50801366-0.89331475j  0.2608311 +0.96414954j  0.2608311 -0.96414954j]

Módulo de las raíces: [1.03809773 1.0098988  1.0098988  1.02454704 1.02454704 1.00138771
 1.00138771 1.0349702  1.0349702  1.01517895 1.01517895 1.03151321
 1.03151321 1.02333621 1.00554483 1.00554483 1.03728888 1.03728888
 1.01332793 1.01332793 1.02766197 1.02766197 0.99880789 0.99880789]

¿Las raíces están fuera del círculo unitario?  False

El modelo no es invertible

también este, pero no es invertible.

In [76]:
modelo=SARIMAX(zpre,
               order=([12],0, [5]),
               seasonal_order=(0,1,[1,2],12)).fit()

modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
Out[76]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX([12], 0, [5])x(0, 1, [1, 2], 12) Log Likelihood -708.534
Date: Sun, 27 Apr 2025 AIC 1427.068
Time: 04:51:22 BIC 1449.437
Sample: 0 HQIC 1435.745
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L12 -0.8858 0.147 -6.011 0.000 -1.175 -0.597
ma.L5 -0.1707 0.044 -3.861 0.000 -0.257 -0.084
ma.S.L12 0.2236 0.140 1.599 0.110 -0.050 0.498
ma.S.L24 -0.6343 0.088 -7.208 0.000 -0.807 -0.462
sigma2 0.5151 0.019 27.358 0.000 0.478 0.552
Ljung-Box (L1) (Q): 10.39 Jarque-Bera (JB): 283.30
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.83 Skew: 0.66
Prob(H) (two-sided): 0.18 Kurtosis: 5.96


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
In [77]:
invertible(modelo)
El polinomio de la parte media móvil es: {'p5': -0.1706643002131486, 'p12': 0.22358740803256344, 'p24': -0.6342671085676491}

El grado del polinomio es: 24

Raíces del polinomio característico: [-1.03809804+0.j         -0.97339404+0.26906235j -0.97339404-0.26906235j
 -0.88525946+0.5157653j  -0.88525946-0.5157653j  -0.71259631+0.70354748j
 -0.71259631-0.70354748j -0.52270685+0.89327583j -0.52270685-0.89327583j
 -0.26081426+0.98110091j -0.26081426-0.98110091j  0.00733501+1.03148763j
  0.00733501-1.03148763j  1.023337  +0.j          0.97341054+0.25216538j
  0.97341054-0.25216538j  0.89999765+0.51572583j  0.89999765-0.51572583j
  0.7125636 +0.72047304j  0.7125636 -0.72047304j  0.50801417+0.89331519j
  0.50801417-0.89331519j  0.26083046+0.96414729j  0.26083046-0.96414729j]

Módulo de las raíces: [1.03809804 1.00989628 1.00989628 1.02454778 1.02454778 1.00138532
 1.00138532 1.0349706  1.0349706  1.01517637 1.01517637 1.03151371
 1.03151371 1.023337   1.00554237 1.00554237 1.03728921 1.03728921
 1.01332536 1.01332536 1.0276626  1.0276626  0.99880555 0.99880555]

¿Las raíces están fuera del círculo unitario?  False

El modelo no es invertible

MODELO PARA LA SERIE TRANSFORMADA¶

In [ ]:
fac = FAC(len(dypre), acf(dypre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos:
r1: 0.1881206416574384
r2: 0.12601374062799692
r6: -0.10023577930193152
r12: -0.40245032162240263
r33: 0.09802046981427195
r45: -0.11985655705785987
In [ ]:
facp = FACP(len(dypre), pacf(dypre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos:
rho 1: 0.18841139999075743
rho 2: 0.09425081653205467
rho 6: -0.0848567418000625
rho 12: -0.4319892618507997
rho 13: 0.08136260875048439
rho 17: -0.11385209887045368
rho 22: -0.08401445717662205
rho 24: -0.2809725408252081
rho 25: 0.10755279636113212
rho 30: -0.09673228306484968
rho 33: 0.1105360450418002
rho 36: -0.14602967562386657
rho 42: -0.11769525648448687
rho 46: -0.09378255882994088
rho 48: -0.16710968828081296
rho 50: 0.09151867173223789
rho 55: 0.09278816335947565
rho 60: -0.09971581325023068
rho 72: -0.13002571917786468
rho 78: -0.07964448430687596
rho 84: -0.07986957194068202
rho 92: -0.13966165907823416
rho 93: 0.11600081821072342
rho 94: -0.1221857229746919
rho 113: -0.08611614522897519
rho 118: -0.08199830520693245
rho 120: -0.11627799482364788

1 MODELO¶

Por la FAC y FACP se cree que puede ser el siguiente modelo:

In [ ]:
modelo=SARIMAX(ypre,
               order=([1, 2, 6],0,1),
               seasonal_order=(1,1,6,12)).fit()

modelo.summary()
Out[ ]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX([1, 2, 6], 0, 1)x(1, 1, [1, 2, 3, 4, 5, 6], 12) Log Likelihood -617.741
Date: Sun, 27 Apr 2025 AIC 1259.482
Time: 18:48:15 BIC 1313.169
Sample: 0 HQIC 1280.309
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.3566 0.144 2.482 0.013 0.075 0.638
ar.L2 0.0211 0.062 0.342 0.732 -0.100 0.142
ar.L6 -0.2754 0.036 -7.653 0.000 -0.346 -0.205
ma.L1 -0.1092 0.149 -0.731 0.465 -0.402 0.184
ar.S.L12 -0.8536 0.246 -3.468 0.001 -1.336 -0.371
ma.S.L12 0.1009 0.245 0.412 0.680 -0.379 0.581
ma.S.L24 -0.7034 0.191 -3.691 0.000 -1.077 -0.330
ma.S.L36 0.0214 0.049 0.435 0.664 -0.075 0.118
ma.S.L48 0.0154 0.054 0.286 0.775 -0.090 0.121
ma.S.L60 0.0169 0.053 0.319 0.750 -0.087 0.121
ma.S.L72 0.0063 0.055 0.114 0.909 -0.101 0.114
sigma2 0.3874 0.020 19.120 0.000 0.348 0.427
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 10.86
Prob(Q): 0.88 Prob(JB): 0.00
Heteroskedasticity (H): 0.83 Skew: -0.20
Prob(H) (two-sided): 0.17 Kurtosis: 3.50


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

RESIDUOS INDEPENDIENTES¶

In [97]:
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
Out[97]:
lb_stat lb_pvalue
1 0.066845 NaN
2 1.537653 NaN
3 1.886894 NaN
4 4.744317 NaN
5 19.994274 NaN
6 20.002045 NaN
7 20.551066 NaN
8 20.563646 NaN
9 25.131429 NaN
10 25.236634 NaN
11 31.550176 NaN
12 31.652621 1.843656e-08
13 31.779184 1.256718e-07
14 31.817201 5.718884e-07
15 31.850384 2.052629e-06
16 34.932704 1.551903e-06
17 42.942269 1.197523e-07
18 48.759895 2.528439e-08
19 50.340543 3.514875e-08
20 50.340543 9.294961e-08
21 51.913020 1.184210e-07
22 52.744140 1.998096e-07
23 57.847442 5.554601e-08
24 57.847451 1.270340e-07
25 59.668596 1.341066e-07
26 59.677107 2.865712e-07
27 59.912799 5.413768e-07
28 60.290463 9.410050e-07
29 60.475407 1.713942e-06
30 69.299462 1.202209e-07
31 69.362452 2.313866e-07
32 69.561924 4.130052e-07
33 73.056696 2.166748e-07
34 74.446295 2.456685e-07
35 76.961663 1.840763e-07
36 76.961687 3.388905e-07
37 78.817354 3.185844e-07
38 80.532525 3.135244e-07
39 83.152029 2.243965e-07
40 83.267861 3.812489e-07
41 84.188435 4.853694e-07
42 86.097859 4.396425e-07
43 86.761902 6.023649e-07
44 86.884845 9.763344e-07
45 90.327760 5.298933e-07
46 91.250891 6.557028e-07
47 98.721616 9.467501e-08
48 98.935868 1.493370e-07

Hay dependencia en los residuos.

RESIDUOS CON MEDIA CERO¶

In [98]:
modelo.resid.mean()
Out[98]:
-0.017262744929477162
In [100]:
ttest_1samp(modelo.resid, 0)
Out[100]:
TtestResult(statistic=-0.6945529062973945, pvalue=0.4875803217531771, df=659)

La media de los residuos puede ser 0.

GRÁFICO DE RESIDUOS¶

In [99]:
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    # line.set_markerfacecolor('green')
    # line.set_markeredgecolor('white')
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

2 MODELO¶

Pero no captura completamente la dependencia temporal, por lo que se agregó la parte completa para la parte AR

In [101]:
modelo=SARIMAX(ypre,
               order=(6,0,1),
               seasonal_order=(1,1,1,12)).fit()

modelo.summary()
Out[101]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(6, 0, 1)x(1, 1, 1, 12) Log Likelihood -602.344
Date: Sun, 27 Apr 2025 AIC 1224.688
Time: 18:50:32 BIC 1269.426
Sample: 0 HQIC 1242.043
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.4302 0.148 2.908 0.004 0.140 0.720
ar.L2 0.0474 0.067 0.702 0.483 -0.085 0.180
ar.L3 -0.0470 0.047 -0.995 0.320 -0.140 0.046
ar.L4 -0.0870 0.043 -2.013 0.044 -0.172 -0.002
ar.L5 -0.1800 0.045 -3.968 0.000 -0.269 -0.091
ar.L6 -0.2197 0.056 -3.904 0.000 -0.330 -0.109
ma.L1 -0.1680 0.154 -1.089 0.276 -0.470 0.134
ar.S.L12 0.0208 0.046 0.451 0.652 -0.070 0.111
ma.S.L12 -0.8907 0.029 -30.470 0.000 -0.948 -0.833
sigma2 0.3657 0.019 19.149 0.000 0.328 0.403
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 5.80
Prob(Q): 0.89 Prob(JB): 0.05
Heteroskedasticity (H): 0.81 Skew: -0.17
Prob(H) (two-sided): 0.12 Kurtosis: 3.32


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

RESIDUOS INDEPENDIENTES¶

In [102]:
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
Out[102]:
lb_stat lb_pvalue
1 0.078029 NaN
2 0.081745 NaN
3 0.110251 NaN
4 0.667736 NaN
5 1.809556 NaN
6 3.209615 NaN
7 3.289842 NaN
8 3.816187 NaN
9 6.473299 NaN
10 7.880155 0.004998
11 9.034590 0.010919
12 9.198592 0.026764
13 9.646686 0.046819
14 9.995740 0.075356
15 10.036647 0.123117
16 10.813321 0.146972
17 20.147979 0.009790
18 24.051508 0.004221
19 24.928482 0.005483
20 25.176289 0.008590
21 25.707594 0.011804
22 30.300176 0.004262
23 32.585774 0.003304
24 32.718036 0.005134
25 33.218651 0.006912
26 33.280245 0.010384
27 33.486398 0.014565
28 33.781554 0.019498
29 33.785406 0.027614
30 39.683798 0.008123
31 39.694743 0.011734
32 39.795172 0.016202
33 43.588079 0.008516
34 43.726951 0.011633
35 45.030235 0.011682
36 45.950944 0.012878
37 49.181595 0.007969
38 50.128815 0.008747
39 52.006371 0.007605
40 52.303949 0.009733
41 53.661655 0.009590
42 55.762429 0.007912
43 56.156343 0.009780
44 56.305459 0.012669
45 60.129401 0.007053
46 63.211636 0.004624
47 69.223074 0.001459
48 69.868377 0.001735

Y el modelo captura correctamente la dependencia temporal.

In [105]:
import modelo_admisible as modadm
In [106]:
modadm.estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.4302180830507938, 'p2': -0.04735397876460526, 'p3': 0.04698747454292404, 'p4': 0.08700716618181399, 'p5': 0.1799715794603173, 'p6': 0.21973455421960683, 'p12': -0.02082020047854175}

Raíces del polinomio característico: [-1.51585027+0.j         -1.16433253+0.66675602j -1.16433253-0.66675602j
 -0.70531602+1.36923083j -0.70531602-1.36923083j  0.00589534+1.27354789j
  0.00589534-1.27354789j  0.89804478+1.3208691j   0.89804478-1.3208691j
  1.62535851+0.j          0.91095431+0.5227982j   0.91095431-0.5227982j ]

Módulo de las raíces: [1.51585027 1.34172793 1.34172793 1.54021549 1.54021549 1.27356154
 1.27356154 1.59724125 1.59724125 1.62535851 1.0503122  1.0503122 ]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)
In [107]:
modadm.invertible(modelo)
El polinomio de la parte media móvil es: {'p1': -0.16798341167047423, 'p12': -0.8906607240754797}

Raíces del polinomio característico: [-1.02312938+0.j         -0.88152548+0.51636208j -0.88152548-0.51636208j
 -0.49857917+0.88672189j -0.49857917-0.88672189j  0.01419601+1.01059688j
  0.01419601-1.01059688j  0.51292531+0.86213479j  0.51292531-0.86213479j
  0.99443323+0.j          0.86733141+0.49151162j  0.86733141-0.49151162j]

Módulo de las raíces: [1.02312938 1.02162467 1.02162467 1.01727917 1.01727917 1.01069658
 1.01069658 1.00317933 1.00317933 0.99443323 0.99691898 0.99691898]

¿Las raíces están fuera del círculo unitario?  False

El modelo no es invertible

3 MODELO¶

In [116]:
modelo=SARIMAX(ypre,
               order=(6,0,0),
               seasonal_order=(1,1,1,12)).fit()

modelo.summary()
Out[116]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(6, 0, 0)x(1, 1, [1], 12) Log Likelihood -603.289
Date: Sun, 27 Apr 2025 AIC 1224.578
Time: 19:07:19 BIC 1264.843
Sample: 0 HQIC 1240.198
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.2721 0.034 7.929 0.000 0.205 0.339
ar.L2 0.1022 0.041 2.473 0.013 0.021 0.183
ar.L3 -0.0255 0.040 -0.644 0.520 -0.103 0.052
ar.L4 -0.0900 0.041 -2.215 0.027 -0.170 -0.010
ar.L5 -0.1974 0.039 -5.011 0.000 -0.275 -0.120
ar.L6 -0.2604 0.038 -6.796 0.000 -0.336 -0.185
ar.S.L12 0.0196 0.046 0.426 0.670 -0.071 0.110
ma.S.L12 -0.8851 0.029 -30.327 0.000 -0.942 -0.828
sigma2 0.3670 0.019 19.130 0.000 0.329 0.405
Ljung-Box (L1) (Q): 0.22 Jarque-Bera (JB): 6.43
Prob(Q): 0.64 Prob(JB): 0.04
Heteroskedasticity (H): 0.80 Skew: -0.19
Prob(H) (two-sided): 0.10 Kurtosis: 3.30


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

PRINCIPIO DE PARSIMONIA¶

In [117]:
parsimonia(modelo) 
Hay coeficientes no significativos, no se cumple el principio de parsimonia

MODELO ADMISIBLE¶

In [118]:
modadm.estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.2721484850682655, 'p2': -0.10222956018637822, 'p3': 0.02552839449395747, 'p4': 0.09002860439631356, 'p5': 0.1974326132169234, 'p6': 0.2604147375399908, 'p12': -0.01958784345380709}

Raíces del polinomio característico: [-1.5324563 +0.j         -0.72632726+1.38599378j -0.72632726-1.38599378j
 -1.13728973+0.63130766j -1.13728973-0.63130766j  0.91165759+1.36955392j
  0.91165759-1.36955392j -0.03290209+1.25835037j -0.03290209-1.25835037j
  1.6720367 +0.j          0.91507128+0.53286918j  0.91507128-0.53286918j]

Módulo de las raíces: [1.5324563  1.56477795 1.56477795 1.30076027 1.30076027 1.64523479
 1.64523479 1.25878045 1.25878045 1.6720367  1.0589169  1.0589169 ]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)
In [119]:
modadm.invertible(modelo)
El polinomio de la parte media móvil es: {'p12': -0.8851100185522583}

Raíces del polinomio característico: [-1.01022217e+00+0.j         -8.74878063e-01+0.50511109j
 -8.74878063e-01-0.50511109j -5.05111085e-01+0.87487806j
 -5.05111085e-01-0.87487806j -5.55111512e-17+1.01022217j
 -5.55111512e-17-1.01022217j  5.05111085e-01+0.87487806j
  5.05111085e-01-0.87487806j  1.01022217e+00+0.j
  8.74878063e-01+0.50511109j  8.74878063e-01-0.50511109j]

Módulo de las raíces: [1.01022217 1.01022217 1.01022217 1.01022217 1.01022217 1.01022217
 1.01022217 1.01022217 1.01022217 1.01022217 1.01022217 1.01022217]

¿Las raíces están fuera del círculo unitario?  True

El modelo es invertible

:)

RESIDUOS INDEPENDIENTES¶

In [120]:
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
Out[120]:
lb_stat lb_pvalue
1 0.001254 NaN
2 0.127263 NaN
3 0.145069 NaN
4 0.355249 NaN
5 1.452317 NaN
6 3.552256 NaN
7 3.845680 NaN
8 4.082125 NaN
9 6.692169 0.009684
10 7.832130 0.019919
11 9.102569 0.027958
12 9.256310 0.055003
13 9.544643 0.089215
14 9.734154 0.136303
15 9.887684 0.195028
16 10.818577 0.212192
17 19.823382 0.019034
18 23.899218 0.007872
19 24.703424 0.010073
20 25.013777 0.014758
21 25.631668 0.019043
22 29.998524 0.007635
23 32.696622 0.005169
24 32.798350 0.007858
25 33.458538 0.009854
26 33.506221 0.014484
27 33.689446 0.019989
28 33.907394 0.026758
29 33.908062 0.037069
30 40.197221 0.010252
31 40.244748 0.014426
32 40.383594 0.019444
33 44.092471 0.010590
34 44.309675 0.013997
35 45.634276 0.013925
36 46.439778 0.015687
37 49.474952 0.010282
38 50.724475 0.010416
39 52.415648 0.009474
40 52.583380 0.012372
41 53.835941 0.012450
42 55.965773 0.010224
43 56.386506 0.012439
44 56.652499 0.015534
45 60.570646 0.008572
46 63.146117 0.006367
47 69.344542 0.001975
48 69.962021 0.002346

Hay dependencia en los primeros dos residuos.

RESIDUOS CON MEDIA CERO¶

In [121]:
modelo.resid.mean()
Out[121]:
-0.03293748637459734
In [122]:
ttest_1samp(modelo.resid, 0)
Out[122]:
TtestResult(statistic=-1.3520655710406162, pvalue=0.1768181651558442, df=659)

La media de los residuos puede ser 0.

RESIDUOS CON VARIANZA CONSTANTE¶

In [123]:
het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))
Out[123]:
(8.59265039527976,
 0.0033752241057990228,
 8.679613399395905,
 0.003331636731645561)

Los residuos tienen varianza constante.

RESIDUOS CON DISTRIBUCIÓN NORMAL¶

In [124]:
print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))
SignificanceResult(statistic=6.693009944604184, pvalue=0.03520718942157648)
(0.02689680139562356, 0.374813270764381)

Los residuos puede que sigan una distribución normal.

In [125]:
normal_relajacion(modelo.resid)
70.00% de los residuos están dentro de ±1σ (esperado ≈ 68%)
95.30% de los residuos están dentro de ±2σ (esperado ≈ 95%)
99.55% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)

Pero son bastante parecidos.

GRÁFICO DE RESIDUOS¶

In [126]:
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    # line.set_markerfacecolor('green')
    # line.set_markeredgecolor('white')
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

MODELO SELECCIONADO¶

Para la serie sin transformar, este es el modelo seleccionado:

$$\text{ARIMA}(13, 0, 0)(0, 1, 4)_{12}$$

Nos quedamos con este modelo porque es el que tiene los p-valores más altos en la prueba de independencia de los residuos. Cumple todos los supuestos, menos el principio de parsimonia, y sus residuos siguen una distribución normal relajada. Tiene un $\text{AIC} = 1367.369$

Para la serie transformada, este es el modelo seleccionado:

$$\text{ARIMA}(6, 0, 0)(1, 1, 1)_{12}$$

Cumple todos los supuestos, menos el principio de parsimonia, y tiene un $\text{AIC} = 1224.578$

In [139]:
modelo = SARIMAX(zpre,
                  order=(13, 0, 0),
                  seasonal_order=(0, 1, 4, 12)).fit()

modelo.summary()
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
Out[139]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(13, 0, 0)x(0, 1, [1, 2, 3, 4], 12) Log Likelihood -665.247
Date: Sun, 27 Apr 2025 AIC 1366.493
Time: 19:21:59 BIC 1447.023
Sample: 0 HQIC 1397.733
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.1420 0.035 4.086 0.000 0.074 0.210
ar.L2 -0.0189 0.037 -0.516 0.606 -0.090 0.053
ar.L3 0.0496 0.037 1.349 0.177 -0.022 0.122
ar.L4 -0.0694 0.042 -1.646 0.100 -0.152 0.013
ar.L5 -0.1131 0.050 -2.274 0.023 -0.211 -0.016
ar.L6 -0.1124 0.060 -1.872 0.061 -0.230 0.005
ar.L7 -0.0152 0.051 -0.299 0.765 -0.115 0.085
ar.L8 -0.0505 0.044 -1.154 0.249 -0.136 0.035
ar.L9 0.0363 0.035 1.030 0.303 -0.033 0.105
ar.L10 -0.0074 0.033 -0.225 0.822 -0.072 0.057
ar.L11 0.0991 0.037 2.677 0.007 0.027 0.172
ar.L12 0.5353 0.082 6.554 0.000 0.375 0.695
ar.L13 -0.0519 0.042 -1.247 0.213 -0.133 0.030
ma.S.L12 -1.3733 0.097 -14.094 0.000 -1.564 -1.182
ma.S.L24 0.3639 0.086 4.217 0.000 0.195 0.533
ma.S.L36 0.0957 0.060 1.599 0.110 -0.022 0.213
ma.S.L48 -0.0740 0.039 -1.888 0.059 -0.151 0.003
sigma2 0.4357 0.028 15.586 0.000 0.381 0.491
Ljung-Box (L1) (Q): 0.11 Jarque-Bera (JB): 344.46
Prob(Q): 0.74 Prob(JB): 0.00
Heteroskedasticity (H): 0.83 Skew: 0.87
Prob(H) (two-sided): 0.18 Kurtosis: 6.12


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
In [138]:
modelo2 = SARIMAX(ypre,
                  order=(6, 0, 0),
                  seasonal_order=(1, 1, 1, 12)).fit()

modelo2.summary()
Out[138]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(6, 0, 0)x(1, 1, [1], 12) Log Likelihood -603.289
Date: Sun, 27 Apr 2025 AIC 1224.578
Time: 19:20:28 BIC 1264.843
Sample: 0 HQIC 1240.198
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.2721 0.034 7.929 0.000 0.205 0.339
ar.L2 0.1022 0.041 2.473 0.013 0.021 0.183
ar.L3 -0.0255 0.040 -0.644 0.520 -0.103 0.052
ar.L4 -0.0900 0.041 -2.215 0.027 -0.170 -0.010
ar.L5 -0.1974 0.039 -5.011 0.000 -0.275 -0.120
ar.L6 -0.2604 0.038 -6.796 0.000 -0.336 -0.185
ar.S.L12 0.0196 0.046 0.426 0.670 -0.071 0.110
ma.S.L12 -0.8851 0.029 -30.327 0.000 -0.942 -0.828
sigma2 0.3670 0.019 19.130 0.000 0.329 0.405
Ljung-Box (L1) (Q): 0.22 Jarque-Bera (JB): 6.43
Prob(Q): 0.64 Prob(JB): 0.04
Heteroskedasticity (H): 0.80 Skew: -0.19
Prob(H) (two-sided): 0.10 Kurtosis: 3.30


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

MÉTODO FORECAST¶

In [ ]:
# Pronóstico para 12 meses hacia el futuro
forecast = modelo.forecast(steps=12)
forecast2 = modelo2.forecast(steps=12)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
  return get_prediction_index(
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception.
  return get_prediction_index(
In [140]:
# Obtener pronóstico con intervalos
pred = modelo.get_forecast(steps=12)
pred2 = modelo2.get_forecast(steps=12)
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
  return get_prediction_index(
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception.
  return get_prediction_index(
In [144]:
# Media pronosticada
forecast_mean = pred.predicted_mean
forecast_mean
Out[144]:
660   -0.276578
661   -0.661192
662   -0.642475
663   -0.132995
664    0.046452
665    0.706953
666    0.536280
667    0.543615
668    0.539949
669   -0.053797
670   -0.541525
671   -0.504177
Name: predicted_mean, dtype: float64
In [ ]:
# Media pronosticada
forecast2_mean = pred2.predicted_mean
forecast2_mean  # El output es un array de numpy, porque tiene la transformacion
Out[ ]:
array([-0.8142261 , -0.58943894, -0.67426909, -0.42378983,  0.14539092,
        0.48758714,  0.50076814,  0.57193519,  0.59488235,  0.21526718,
       -0.35086441, -0.49899387])
In [142]:
# aplicar la inversa de la estandarización
# Calcular media y desviación estándar de la serie original
mu = np.mean(pre)
sigma = np.std(pre, ddof=0)  # ddof=0 porque así lo usa scipy.stats.zscore por default
In [ ]:
# Desestandarizar
forecast_mean = forecast_mean * sigma + mu 
forecast_mean 
Out[ ]:
660     540.502158
661     201.870616
662     218.349640
663     666.919279
664     824.911936
665    1406.446553
666    1256.178392
667    1262.636919
668    1259.408611
669     736.647934
670     307.230864
671     340.113522
Name: predicted_mean, dtype: float64
In [146]:
forecast2_mean = forecast2_mean * sigma + mu
forecast2_mean
Out[146]:
array([  67.13264749,  265.04529754,  190.35704854,  410.89015873,
        912.02229063, 1213.30710951, 1224.91224955, 1287.57089294,
       1307.77459904,  973.54446573,  475.09695852,  344.67717081])
In [84]:
# Intervalos de confianza
int_conf = pred.conf_int()
int_conf
Out[84]:
lower y upper y
660 -1.566687 1.023425
661 -1.969017 0.645214
662 -1.962076 0.652155
663 -1.451584 1.165524
664 -1.273046 1.348105
665 -0.629675 2.012863
666 -0.807778 1.859318
667 -0.787023 1.882975
668 -0.791284 1.883888
669 -1.390333 1.285407
670 -1.872244 0.804378
671 -1.848288 0.850488
In [85]:
# Desestandarizar los intervalos de confianza
int_conf = int_conf * sigma + mu
int_conf
Out[85]:
lower y upper y
660 -595.367749 1685.082779
661 -949.596462 1352.088860
662 -943.485554 1358.199735
663 -494.026019 1810.193098
664 -336.832977 1970.945813
665 229.619785 2556.228270
666 72.809546 2421.040094
667 91.083533 2441.868542
668 87.331818 2442.672581
669 -440.097428 1915.743419
670 -864.393291 1492.223637
671 -843.301917 1532.821294
In [86]:
# Asignar los índices de fecha a los intervalos de confianza
forecast_index = pd.date_range(start=pre.index[-2], periods=13, freq='ME')[1:]
In [87]:
plt.figure()

ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3))  # Pone cada 3 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha

# Serie original
plt.plot(pre_total, label='Observado')

# Pronóstico
plt.plot(forecast_index, forecast_mean, label='Pronóstico', color=sns.color_palette()[1])

# Intervalos
plt.fill_between(forecast_index,
                    int_conf.iloc[:, 0], int_conf.iloc[:, 1], 
                    color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")

plt.legend()
plt.title('Método Forecast para modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
plt.savefig('imagenes/05-01-metodo-forecast.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image
In [88]:
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4))  # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y'))   # Formato de fecha

# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2000-01-01']

# Plotear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, forecast_mean , label='Pronóstico', color=sns.color_palette()[1], linewidth=3, alpha=0.7)

# Intervalos
plt.fill_between(forecast_index,
                    int_conf.iloc[:, 0], int_conf.iloc[:, 1], 
                    color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")


plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Método Forecast para Modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
# plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/05-02-metodo-forecast-zoom.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

PRONÓSTICO ÓPTIMO¶

Se calcula en 04-Calculo-del-Pronostico-Optimo

In [89]:
pre_total.to_csv('04-01-Serie-Original.csv', index=False) # Guardar la serie original
modelo.resid.to_csv('04-02-Serie-Residuos.csv', index=False) # Guardar la serie de residuos
In [90]:
modelo.summary()
Out[90]:
SARIMAX Results
Dep. Variable: y No. Observations: 660
Model: SARIMAX(13, 0, 0)x(0, 1, [1, 2, 3, 4], 12) Log Likelihood -665.685
Date: Sun, 27 Apr 2025 AIC 1367.369
Time: 04:52:23 BIC 1447.899
Sample: 0 HQIC 1398.609
- 660
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 0.1369 0.035 3.963 0.000 0.069 0.205
ar.L2 -0.0190 0.036 -0.522 0.602 -0.090 0.052
ar.L3 0.0501 0.037 1.363 0.173 -0.022 0.122
ar.L4 -0.0697 0.042 -1.662 0.097 -0.152 0.013
ar.L5 -0.1117 0.050 -2.254 0.024 -0.209 -0.015
ar.L6 -0.1101 0.060 -1.839 0.066 -0.227 0.007
ar.L7 -0.0104 0.051 -0.205 0.838 -0.110 0.089
ar.L8 -0.0511 0.044 -1.170 0.242 -0.137 0.034
ar.L9 0.0331 0.035 0.944 0.345 -0.036 0.102
ar.L10 -0.0093 0.033 -0.281 0.778 -0.074 0.055
ar.L11 0.1034 0.037 2.802 0.005 0.031 0.176
ar.L12 0.5385 0.080 6.699 0.000 0.381 0.696
ar.L13 -0.0460 0.041 -1.115 0.265 -0.127 0.035
ma.S.L12 -1.3756 0.100 -13.725 0.000 -1.572 -1.179
ma.S.L24 0.3674 0.086 4.275 0.000 0.199 0.536
ma.S.L36 0.0973 0.060 1.635 0.102 -0.019 0.214
ma.S.L48 -0.0784 0.039 -1.994 0.046 -0.155 -0.001
sigma2 0.4340 0.030 14.503 0.000 0.375 0.493
Ljung-Box (L1) (Q): 0.04 Jarque-Bera (JB): 347.08
Prob(Q): 0.85 Prob(JB): 0.00
Heteroskedasticity (H): 0.83 Skew: 0.87
Prob(H) (two-sided): 0.17 Kurtosis: 6.14


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
In [102]:
modelo.params
Out[102]:
ar.L1       0.136917
ar.L2      -0.019014
ar.L3       0.050067
ar.L4      -0.069698
ar.L5      -0.111714
ar.L6      -0.110100
ar.L7      -0.010426
ar.L8      -0.051137
ar.L9       0.033088
ar.L10     -0.009268
ar.L11      0.103359
ar.L12      0.538506
ar.L13     -0.046024
ma.S.L12   -1.375601
ma.S.L24    0.367441
ma.S.L36    0.097350
ma.S.L48   -0.078355
sigma2      0.434001
dtype: float64
$$ \text{ARIMA}(13,0,0) \times (0,1,4)_{12} \text{ para $X_t$} $$

con coeficientes:

  • $\phi_1 = -0.136917$
  • $\phi_2 = 0.019014$
  • $\phi_3 = -0.050067$
  • $\phi_4 = 0.069698$
  • $\phi_5 = 0.111714$
  • $\phi_6 = 0.110100$
  • $\phi_7 = 0.010426$
  • $\phi_8 = 0.051137$
  • $\phi_9 = -0.033088$
  • $\phi_{10} = 0.009268$
  • $\phi_{11} = -0.103359$
  • $\phi_{12} = -0.538506$
  • $\phi_{13} = 0.046024$

y para la parte estacional MA:

  • $\Theta_1 = 1.375601$
  • $\Theta_2 = -0.367441$
  • $\Theta_3 = -0.097350$
  • $\Theta_4 = 0.078355$
$$ \Phi_0(B^{12}) \, \phi_13(B) \, \nabla^0 \nabla_{12}^1 X_t = \Theta_4(B^{12}) \, \theta_0(B) \, \varepsilon_t, $$

Donde $\varepsilon_t \sim \mathcal{N}(0, 1)$

Desarrollo para el pronóstico óptimo¶

Del lado izquierdo:

$$ (1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_{13} B^{13})(X_t - X_{t-12}) = X_t - X_{t-12} - \phi_1 (X_{t-1} - X_{t-13}) - \phi_2 (X_{t-2} - X_{t-14}) - \cdots - \phi_{13}(X_{t-13} - X_{t-25}) $$$$ = X_t - X_{t-12} - \sum_{h=1}^{13} \phi_k (X_{t-h} - X_{t-h-12}) $$

Del lado derecho:

$$ (1 - \Theta_1 B^{12} - \Theta_2 B^{24} - \Theta_3 B^{36} - \Theta_4 B^{48})\varepsilon_t = \varepsilon_t - \Theta_1 \varepsilon_{t-12} - \Theta_2 \varepsilon_{t-24} - \Theta_3 \varepsilon_{t-36} - \Theta_4 \varepsilon_{t-48} $$

despejando $X_t$:

$$ X_t = X_{t-12} + \sum_{k=1}^{13} \phi_k (X_{t-k} - X_{t-k-12}) + \varepsilon_t - \Theta_1 \varepsilon_{t-12} - \Theta_2 \varepsilon_{t-24} - \Theta_3 \varepsilon_{t-36} - \Theta_4 \varepsilon_{t-48} $$

Por lo que, el pronóstico óptimo sería:

$$ X_t(h) = \underset{t}{\mathbb{E}}[X_{t+h-12}] + \sum_{k=1}^{13} \phi_k \left( \underset{t}{\mathbb{E}}[X_{t+h-k}] - \underset{t}{\mathbb{E}}[X_{t+h-k-12}] \right) + \underset{t}{\mathbb{E}}[\varepsilon_{t+h}] - \sum_{j=1}^4 \Theta_j \, \underset{t}{\mathbb{E}}[\varepsilon_{t+h-12j}] $$
In [91]:
# Intervalo de Predicción
# Pronosntico optimo[i] +- 1.96*(suma de psis)^{1/2}*std(residuos)

# Calcula de las espis
psi0 = -1

# psi[i] = teta[i] + Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d] si i=1,2,...,q
# psi[i] = Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d]           si i>q
In [92]:
# pronostico optimo 04-Calculo-del-Pronostico-Optimo

# 893.2021421
# 218.8353082
# 14.55589463
# 46.06671573
# 102.6854446
# 909.4688428
# 328.2252003
# 431.4633135
# 1179.309604
# 567.0569539
# 460.930121
# 514.9098934

# Guardarlos en una serie

pronostico_optimo = pd.Series([893.2021421, 218.8353082, 14.55589463, 46.06671573, 102.6854446, 909.4688428, 328.2252003, 431.4633135, 1179.309604, 567.0569539, 460.930121, 514.9098934],
                            index=pd.date_range(start='2009-01-01', periods=12, freq='ME'))
In [93]:
# Intervalo de Predicción
# Pronosntico optimo[i] +- 1.96*(suma de psis)^{1/2}*std(residuos)

# Calcula de las espis
psi0 = -1

# psi[i] = teta[i] + Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d] si i=1,2,...,q
# psi[i] = Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d]           si i>q

phi1 = modelo.params[0]
teta12 = modelo.params[1]

psi1 = phi1*psi0
psi2 = phi1*psi1
psi3 = phi1*psi2 
psi4 = phi1*psi3
psi5 = phi1*psi4
psi6 = phi1*psi5
psi7 = phi1*psi6
psi8 = phi1*psi7
psi9 = phi1*psi8
psi10 = phi1*psi9
psi11 = phi1*psi10
psi12 = phi1*psi11 + teta12
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\294999966.py:10: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  phi1 = modelo.params[0]
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\294999966.py:11: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  teta12 = modelo.params[1]
In [94]:
psis = [psi0, psi1, psi2, psi3, psi4, psi5, psi6, psi7, psi8, psi9, psi10, psi11, psi12]
In [95]:
sigmax = np.std(modelo.resid, ddof=0)  # ddof=0 porque así lo usa scipy.stats.zscore por default
In [96]:
# 1.96 es el valor crítico para un intervalo de confianza del 95%

# primer intervalo
u0 = pronostico_optimo[0] + 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l0 = pronostico_optimo[0] - 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# segundo intervalo
u1 = pronostico_optimo[1] + 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l1 = pronostico_optimo[1] - 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# tercer intervalo
u2 = pronostico_optimo[2] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l2 = pronostico_optimo[2] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# cuarto intervalo
u3 = pronostico_optimo[3] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l3 = pronostico_optimo[3] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# quinto intervalo
u4 = pronostico_optimo[4] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l4 = pronostico_optimo[4] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# sexto intervalo
u5 = pronostico_optimo[5] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l5 = pronostico_optimo[5] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# septimo intervalo
u6 = pronostico_optimo[6] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l6 = pronostico_optimo[6] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# octavo intervalo
u7 = pronostico_optimo[7] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l7 = pronostico_optimo[7] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# noveno intervalo
u8 = pronostico_optimo[8] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l8 = pronostico_optimo[8] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# decimo intervalo
u9 = pronostico_optimo[9] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l9 = pronostico_optimo[9] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# decimo primero intervalo
u10 = pronostico_optimo[10] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l10 = pronostico_optimo[10] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# decimo segundo intervalo
u11 = pronostico_optimo[11] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l11 = pronostico_optimo[11] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:4: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u0 = pronostico_optimo[0] + 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:5: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l0 = pronostico_optimo[0] - 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:8: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u1 = pronostico_optimo[1] + 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:9: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l1 = pronostico_optimo[1] - 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:12: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u2 = pronostico_optimo[2] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:13: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l2 = pronostico_optimo[2] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:16: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u3 = pronostico_optimo[3] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:17: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l3 = pronostico_optimo[3] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:20: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u4 = pronostico_optimo[4] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:21: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l4 = pronostico_optimo[4] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:24: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u5 = pronostico_optimo[5] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:25: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l5 = pronostico_optimo[5] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:28: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u6 = pronostico_optimo[6] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:29: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l6 = pronostico_optimo[6] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:32: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u7 = pronostico_optimo[7] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:33: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l7 = pronostico_optimo[7] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:36: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u8 = pronostico_optimo[8] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:37: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l8 = pronostico_optimo[8] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:40: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u9 = pronostico_optimo[9] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:41: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l9 = pronostico_optimo[9] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:44: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u10 = pronostico_optimo[10] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:45: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l10 = pronostico_optimo[10] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u11 = pronostico_optimo[11] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:49: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l11 = pronostico_optimo[11] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
In [97]:
# Guardar los resultados en un dataframe
pronostico_optimo_df = pd.DataFrame({
    'Pronóstico': pronostico_optimo,
    'Límite Inferior': [l0, l1, l2, l3, l4, l5, l6, l7, l8, l9, l10, l11],
    'Límite Superior': [u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11]
}, index=pronostico_optimo.index)
In [98]:
import matplotlib.dates as mdates

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4))  # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y'))   # Formato de fecha

# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2000-01-01']

# Plottear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, pronostico_optimo, label='Pronóstico', color=sns.color_palette()[1], linewidth=3, alpha=0.7)

# Intervalos
plt.fill_between(forecast_index,
                    pronostico_optimo_df['Límite Inferior'], pronostico_optimo_df['Límite Superior'], 
                    color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")

plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Pronóstico Óptimo con Esperanzas Condicionales')
plt.ylim(-1000, 6000)
plt.grid(True)
plt.savefig('imagenes/05-03-pronostico-optimo.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

COMPARACIÓN¶

In [107]:
# Comparar el pronóstico óptimo con el pronóstico del modelo SARIMAX

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4))  # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y'))   # Formato de fecha

# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2005-01-01']

# Plottear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, pronostico_optimo, label='Pronóstico Óptimo', color=sns.color_palette()[1], linewidth=3, alpha=0.7)
plt.plot(forecast_index, forecast_mean, label='Método Forecast', color=sns.color_palette()[2], linewidth=3, alpha=0.7)

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")

plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Pronóstico Óptimo vs Método Forecast')
plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/05-04-comparacion.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image

ACTUALIZACIÓN DE PRONÓSTICO¶

In [100]:
# Conseguimos la primera observación de la serie test
valornuevo = pre_test[0]

# Estandarizamos
znv = (valornuevo - mu) / sigma

# Crear una nueva marca de tiempo para la nueva observación (usando el siguiente período después de la última fecha de la serie estandarizada)
new_date = zpre.index[-1] + pd.DateOffset(months=1)
new_obs = pd.Series(znv, index=[new_date])

# Actualizar la serie estandarizada con la nueva observación (si es necesario) utilizando pd.concat
zpre_updated = pd.concat([zpre, new_obs])

# Actualizar el modelo con la nueva observación sin re-calibrarlo
modelo_updated = modelo.append([new_obs], refit=False)

# Generar pronóstico para los próximos 12 meses (en escala estandarizada)
forecast_std_updated = modelo_updated.forecast(steps=11)

# Transformar el pronóstico de vuelta a la escala original
forecast_updated = forecast_std_updated * sigma + mu

# Asignar los índices de fecha a los intervalos de confianza
forecast_index2 = pd.date_range(start=pre.index[-1], periods=12, freq='ME')[1:]
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\2112643297.py:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  valornuevo = pre_test[0]
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
  return get_prediction_index(
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception.
  return get_prediction_index(
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
  return get_prediction_index(
In [101]:
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4))  # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y'))   # Formato de fecha

# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2005-01-01']

# Plotear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index2, forecast_updated , label='Pronóstico Actualizado', color=sns.color_palette()[1], linewidth=3, alpha=0.7)
plt.plot(forecast_index, forecast_mean , label='Pronóstico Original', color=sns.color_palette()[2], linewidth=3, alpha=0.7)

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")


plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Método Forecast para Modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
# plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/06-01-Actualizacion-de-pronostico-metodo-forecast.svg', bbox_inches='tight')
plt.show()
No description has been provided for this image